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ABSTRACT 

p j The formation and evolution of protoplanetary systems, the breeding grounds of planet 

rj^j formation, is a complex dynamical problem that involves many orders of magnitudes. 

To serve this purpose, we present a new hybrid algorithm that combines a Fokker- 
i-£h Planck approach with the advantages of a pure direct-summation TV— body scheme, 

with a very accurate integration of close encounters for the orbital evolution of the 
Q larger bodies with a statistical model, envisaged to simulate the very large number of 

S— ( smaller planetesimals in the disc. Direct-summation techniques have been historically 

developped for the study of dense stellar systems such as open and globular clus- 
ters and, within some limits imposed by the number of stars, of galactic nuclei. The 
number of modifications to adapt direct-summation N— body techniques to planetary 
dynamics is not undemanding and requires modifications. These include the way close 
J> encounters are treated, as well as the selection process for the "neighbour radius" of the 

■^J- particles and the extended Hermite scheme, used for the very first time in this work, 

Cf\ as well as the implementation of a central potential, drag forces and the adjustment of 

the regularisation treatment. For the statistical description of the planetesimal disc we 
^*>0 employ a Fokker-Planck approach. We include dynamical friction, high- and low-speed 

l/-} encounters, the role of distant encounters as well as gas and collisional damping and 

then generalise the model to inhomogenous discs. We then describe the combination 
t-H of the two techniques to address the whole problem of planetesimal dynamics in a 

t-H realistic way via a transition mass to integrate the evolution of the particles according 

^■J to their masses. 
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1 INTRODUCTION 

The formation of a planetary system is closely related to 
the formation of the host star itself. Cool molecular clouds 
collapse and fragment into smaller substructures which are 
the seeds for subsequent star formation. Angular momen- 
tum conservation in the forming clumps forces the infalling 
matter into a disc-like structure. The subsequent viscous 
evolution of the disc leads to a transport of angular momen- 
tum which channels gas to the protostar in the centre. These 
protoplanetary discs are the birth place of planets (for a de- 
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tailed review see Armitage 2010). Embedded dust grains are 
the seed for the enormous growth to bodies of planetary size. 
The first hint to the structure of protoplanetary discs has 
been provided by our own solar system. Through "smearing 
out" all planets and adding the missing fraction of volatile 
elements, one can estimate the structure and mass of the 
protoplanetary disc. Since the efficiency of planet formation 
is unknown, this yields only a lower limit - the minimum 
mass solar nebula (Hayashi 1981). The inferred surface den- 
sity decreases with the distance from the sun as oc r~ 3//2 . 

Further insight has been obtained by the detection of an 
infrared excess of some stars, i. e. additional infrared radia- 
tion that does not originate from the star but an unresolved 
disc. Advancements in observation led in the mid-90s to the 
direct imaging of nearby star forming regions which opened a 
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new nourishing field in astronomy (see O'dell et al. 1993, for 
an example with Hubble Space Telescope images)). Since 
then a big amount of observations of protoplanetary discs 
across the electromagnetic range has been obtained and in- 
terpreted, both space-based (ISO, Spitzer and now Herschel) 
and in the ground (VLT, VLT-I, Subaru, Keck). We refer the 
reader to "The Star Formation Newsletter" URL 1 for an 
overview of the recent papers in star- and planet formation 
and the book review "Protostars and Planets V" (Reipurth 
et al. 2007). 

Protoplanetary discs masses cover a range from 10~ 3 
to O.1M with a peak around O.O1M (data from Tau- 
rus/Ophiuchus Beckwith 1996), in accordance with mass es- 
timates deduced from the minimum mass solar nebula. Since 
the disc lifetime can not be measured directly, it is derived 
indirectly from the age of young, naked (i. e. discless) stars 
which sets an upper limit. Pre-main sequence evolutionary 
tracks are used to gauge the stellar ages, providing lifetimes 
of a few 10 6 years. The subsequent evolution of the disc 
proceeds in several stages. 

Two different scenarios have been proposed to explain 
the formation of kilometre-sized planetesimals. 

(i) One process is the gravitational instability of the dust 
component in a protoplanetary disc that leads to the direct 
formation of larger bodies. Goldreich &: Ward (1973) pro- 
pose that an initial growth phase of dust grains leads to a 
thin dust disc that undergoes a gravitational collapse. As 
the dense dust layer decouples from the gas, it rotates with 
the local Keplerian velocity, whereas the gas component ro- 
tates slower as it is partially pressure supported. This gives 
rise to a velocity shear at the boundary, which may excite 
turbulence through the Kelvin-Helmholtz-instability. Since 
the motion of small dust grains in the boundary layer is cou- 
pled to the gas, the turbulent velocity field could suppress 
the formation of a stratified dust layer, which is a necessary 
prerequisite for the gravitational instability (Weidenschilling 
1977). 

(ii) The collisional agglomeration of dust particles is an 
opposed formation mechanism. Relative velocities are dom- 
inated by the Brownian motion in the early phases of the 
growth process. This mechanism becomes increasingly in- 
efficient with growing mass, but successively the particles 
decouple from the gas and settle to the midplane - a pro- 
cess that yields even larger velocities with increasing mass. 
The sedimentation initiates a growth mode that is similar 
to the processes in rain clouds: Larger grains drop faster, 
thus accreting smaller grains on their way to the midplane. 
Turbulence may modify this basic growth scenario by forc- 
ing the dust grains in a convection-like motion. Dust grains 
still grow during the settling process, but the turbulent ve- 
locity field could mix up dust from the midplane, and a new 
cycle begins. Each cycle adds a new layer to the dust grains 
- a mechanism that also operates in hail clouds - until the 
grains are large enough to decouple from the turbulent mo- 
tion. Again, turbulence plays an important role in determin- 
ing the growth mode and the relative velocities. While the 
relative velocities are high enough to allow for a fast growth, 
it is not clear a priory that collisions are sticky enough to 

1 http : //www . if a . hawaii . edu/users/reipurth/newsletter . 
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allow for a net growth. High speed encounters lead to frag- 
mentation, which counteracts agglomeration (e.g. Blum & 
Wurm 2000, and references therein). An important bottle- 
neck in the agglomeration process is the fast orbital decay 
of 1 m sized boulders. Their orbital lifetime is as short as 
100 years, and a quick increase in size - at least over one 
order of magnitude - is needed to reduce the radial drift 
significantly. 

To overcome the difficulties associated with each of these 
scenarios, modifications have been proposed. Magneto- 
hydrodynamic simulations include electro-magnetic interac- 
tions in hydrodynamical calculations. See the reviews of Bal- 
bus & Hawley 1998 and Balbus. 2003. The MHD simulations 
by Johansen et al. (2006) show that trapping of larger parti- 
cles in turbulent vortices helps in increasing the orbital life- 
time, but could also trigger local instabilities that may lead 
to the direct formation of planetesimals (Inaba et al. 2005). 
Johansen et al. (2007) describe a gravoturbulence mecha- 
nism as a feasible pathway to planetesimal formation in ac- 
creting circumstellar discs. 

The details of agglomeration have drawn a lot of at- 
tention and are still under question (see Kempf et al. 1999; 
Paszun & Dominik 2009; Wada et al. 2009), but the succes- 
sive agglomeration of planetesimals is commonly accepted. 

1.1 Formation of protoplanets 

The further growth of planetesimals proceeds through mu- 
tual collisions, where the initial phase involves a large num- 
ber of particles and is well described by a coagulation equa- 
tion (Safronov 1969). While earlier works (e.g. Nakagawa 
et al. 1983) focussed mainly on the evolution of the size dis- 
tribution, subsequent refinements of the evolution of the ran- 
dom velocities showed that it is important to evolve the size 
distribution and the velocity dispersion in a consistent way. 
A fixed velocity dispersion is an oversimplification, which 
changes the growth mode and increases the growth timescale 
as well (Wetherill 1989). 

The initial growth is quite democratic. All planetesimals 
grow roughly at the same rate and the maximum of the size 
distribution is shifted gently towards larger sizes. As soon 
as gravitational focusing and dynamical friction become im- 
portant, the growth mode changes to a qualitatively differ- 
ent mechanism. Efficient gravitational focusing leads to a 
growth timescale (which we denote as M/M) that decreases 
with mass. Hence larger bodies grow faster than the smaller 
planetesimals, a trend that is further supported by energy 
cquipartition due to planetesimal-planetesimal encounters. 
This dynamical friction keeps the largest bodies on nearly 
circular orbit, thus the relative velocities are small and gravi- 
tational focusing remains efficient. Smaller planetesimals are 
stirred up into eccentricity orbits, which slows down their 
growth rate compared to the largest bodies. This acceler- 
ated growth, denoted as runaway growth of planetesimals 
(see e.g Greenberg et al. 1978; Wetherill 1989, 1993), short- 
ens the growth timescale to a few 10 5 years. The term run- 
away growth stresses that the growth timescale of a particle 
decreases with mass, hence the largest body "runs away" to 
the high mass end of the distribution (see Kokubo 1995). 

While energy equipartition increases the velocity dis- 
persion with decreasing mass of the planetesimals, addi- 
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tional damping due to the gaseous disc leads to a turn-over 
at smaller sizes. However, higher relative velocities may lead 
nevertheless to destructive encounters, but these fragmen- 
tation events could even speed up the growth (Chambers 
2006; Bromley & Kenyon 2006, 2010). Since smaller bodies 
are more subjected to gas drag, their velocity dispersion is 
smaller which allows an efficient accretion by the runaway 
bodies. Moreover, smaller particles damp the velocity dis- 
persion of the largest bodies more efficiently. 

As runaway proceeds, the system becomes more and 
more dominated by few big bodies - the protoplanets . Due 
to the dominance of few, very large bodies one can not use 
statistical methods anymore to simulate the problem. 

Kobayashi et al. (2010) derived analytical expressions 
for the final masses of these planetary "embryos", as they 
call them, including the role of planetesimal depletion due to 
collisional disruption. They conclude that the final mass in 
the minimum-mass solar nebula at several AU can achieve 
~ 0.1 Earth mass within 10 7 years. 

1.2 Oligarchic growth 

The runaway growth of large bodies (i. e. protoplanets) 
ceases to be efficient as soon as the protoplanets start to 
control the velocity dispersion of the remaining planetesi- 
mals in their vicinity. Gravitational focusing becomes less 
effective, therefore the growth timescale increases with size 
and the growth mode changes to oligarchic growth. The pro- 
toplanets still grow faster than the field planetesimals 2 , but 
the masses of the protoplanets remain comparable. 

A combination of dynamical friction due to the field 
planetesimals and perturbations from the neighbouring pro- 
toplanets conserves a separation of five to ten Hill radii 
between neighbouring bodies. Therefore only planetesimals 
from a limited area, the feeding zone, are accreted by a given 
protoplanet. If this zone is emptied, they have reached their 
final isolation mass (Kokubo et al. 1998; Kokubo & Ida 2000, 
2002). 

As the damping of the remaining field planetesimals is 
weak enough, further growth is dominated by mutual per- 
turbations among the protoplanets, which leads to giant im- 
pacts. Protoplanets beyond the "snow line" (or ice conden- 
sation point) can grow larger than 5-15 earth masses and 
initiate the formation of giant planets. If the protoplanetary 
disc is very massive in the inner planetary system, this may 
lead to an in-situ formation of hot Jupiters (see Bodenheimer 
et al. 2000). 

Kokubo et al. (2006) used the oligarchic growth model 
of protoplanets to statistically quantify the giant impact 
stage (i.e. collisions of protoplanets to form planets) with 
N— body simulations. They found that for steeper surface 
density profiles, large planets usually form closer to the star. 

1.3 Migration 

The proposed three-stage scenario of planet formation cov- 
ers the dominant growth processes, but a major mechanism 

2 The term "field planetesimals" denotes in the following the 
smooth component of smaller planetesimals. 



is still missing - the migration of bodies in the system. Mi- 
gration is a generic term that summarises a set of different 
mechanisms that lead to secular radial drift of bodies (see 
e.g. the review of Papaloizou 2006; Armitage 2010): 

(i) The dissipation due to the remaining gas disc leads 
to an orbital decay of the planetesimals. While this poses a 
severe problem for 1 m-sized objects, larger bodies drift very 
slowly inward. One denotes this process as type migration. 

(ii) Planets which are embedded in a gaseous disc launch 
spiral density waves at the inner and outer Lindblad reso- 
nances, which leads to an exchange of angular momentum 
with the resonant excited waves. This type I migration leads 
to a robust inward migration independent of the density pro- 
file. The perturbation from the planet is small, hence linear 
perturbation theory is in principle applicable (Ward 1986), 
but recent work proves that non-linear and radiative effects 
can be quite important (Paardekooper et al. 2011). 

(iii) If the protoplanet is massive enough, it opens a gap 
in the gaseous disc and excites waves through tidal interac- 
tion with the gaseous disc. An imbalance of the exchange 
of angular momentum with the inner and outer part of the 
disc leads to type II migration. The strong interaction be- 
tween planet and disc leaves the linear regime and requires 
a numerical solution of the hydrodynamic equations (Lin & 
Papaloizou 1979). 

(iv) Even with an opened gap, the planet still channels 
gas between the inner and outer part of the disc. While this 
corotational flow is to some extent already present during 
type II migration, it dominates the angular momentum bal- 
ance in the case of type III migration due to an asymmetry 
in the leading and trailing part of the flow. An imbalance 
allows for an efficient exchange of angular momentum with 
this corotational gas stream, which gives rise to a remark- 
ably fast migration (Masset & Papaloizou 2003). Type III 
migration is subjected to a positive feedback: A faster mi- 
gration increases the asymmetry in the corotational flow, 
which speeds up the migration. Both inward and outward 
migration are possible. 

The four migration types modify the three-stage scenario in 
different ways. 

Type migration is most efficient for small planetesimals 
(i. e. 1 km or smaller). It poses a severe problem during the 
early stages of planet formation, as it may induce a signifi- 
cant loss of solid material, accompanied by a global change 
in the initial surface density (Kornet et al. 2001). The im- 
portance of this process diminishes as planetesimal growth 
proceeds, but it still leads to the loss of collisional fragments 
during the final disc clearing. 

Type I is most important for protoplanets (i.e. 0.1 earth 
masses or larger). It leads to an orbital decay of protoplan- 
ets, but this does not only imply a loss of protoplanets, but 
also breaks the conditions for isolation. Migrating bodies can 
accrete along their way through the disc and are thus not 
constrained to a fixed feeding zone, which may increase the 
isolation mass. 

Type II and type III mainly influences giant planets (larger 
than 10 earth masses) causing an inward or outward migra- 
tion, depending on the angular momentum exchange bal- 
ance. It can explain the large number of giant planets close 
to their host star (hot Jupiters) found in extrasolar systems. 
An important issue is the timescale of the migration process. 
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If the migration is too fast, virtually all planets spiral inward 
and leave an empty system behind. 

Migration is a powerful process with the capability to 
reshape an entire planetary system, at least for gas giants, 
since terrestrial planets are likely not affected by migration 
in gaseous discs (Bromley & Kenyon 2011). However, it also 
requires some "parking mechanism" which terminates mi- 
gration before all planets (or protoplanets) are lost. Inho- 
mogeneities in the gaseous disc may change the crucial mo- 
mentum balance of the inner and outer part of the disc, thus 
stalling or even reversing the drift of a planet. 

The migration processes end after the dissipation of the 
gaseous disc due to photoevaporation or star-star encounters 
(i. e. after a few 10 6 to 10 7 years). 

The formation of a planetary system is a vital pro- 
cess that is driven by the interplay between the different 
growth phases and the migration of planets and protoplan- 
etary cores (i.e. the precursors of giant planets). While the 
preceding sections only summarised the main evolutionary 
processes, even more processes could influence the forma- 
tion of planetary systems. A fast accretion of giant planets 
in the outer parts of a planetary system could introduce fur- 
ther perturbations on the inner part and may even trigger 
the formation of terrestrial planets. Moreover, the stellar en- 
vironment in dense star clusters and multiple stellar systems 
also perturbs planet formation, which would require an even 
broader view on the problem. Last but not least, the exci- 
tation due to gravitational coupling to gas surface density 
fluctuations plays an important role (see e.g. Ida et al. 2008). 

Any approach to planet formation can hardly include 
this wealth of different phenomena, thus it is important to 
focus on a well-defined subproblem. In this work we focus 
on the formation of protoplanets for the following reasons: 

(i) The size, growth timescale and spacing of the proto- 
planets is a key element in the planet formation process. 

(ii) The protoplanet growth is well-defined by different 
growth modes. It starts with the already formed planetes- 
imals (« 1-5 km) and ends, when orbital crossing of the 
protoplanets initiates the final growth phase. 

(iii) The planetesimals are large enough to treat the re- 
maining gaseous disc as a small perturbation. 

(iv) The protoplanets are small enough to neglect tidal 
interaction with the disc in the inner planetary system. Col- 
lective planetesimal-protoplanet interaction are also negli- 
gible. 

Though the protoplanet formation is a well-posed subprob- 
lem, our approach has to incorporate various mechanisms 
and techniques to cover the full size range of the problem 
(see the review of Armitagc 2010, for a extensive review of the 
problem and references therein). However, it is still accessi- 
ble to theoretical calculations to some extent which provide 
a safe ground for the analysis of the results. 

In next section 2 we briefly summarise the analytical 
model that we use for our initial models. Then we intro- 
duce the problem that we will address later, as well as the 
timescales and quantities of interest for the general prob- 
lem of planetesimal formation: In sections 3, 4, 5 we give a 
description for a test particle moving around a body with 
a central mass M, the Hill's problem and the protoplanet 
growth from a theoretical point of view, respectively. 



Subsequently, in section 6 and its corresponding subsec- 
tions, we ascertain the direct-summation part of the general 
algorithm, as well as the required modifications and addi- 
tions to tailor it to the specific planetesimal problem. In 
section 9 we introduce our collisional treatment, including 
fragmentation, collisional cascades, migration and coagula- 
tion. The statistical part of the algorithm is described in 
detail in section 16. Finally, in section 17 we explain how 
to bring together both schemes into a single numerical tool 
and in section 18 we give our conclusions and compare our 
scheme to other numerical approaches. 



2 INITIAL MODELS 

The basis for all planet formation models is the structure of 
the protoplanetary disc. We summarise the pioneering work 
of Hayashi (1981) to have a robust initial model at hand. 
Subsequent evolution of the disc may change this simple 
approach, but it is still a valuable guideline. 

A basic estimate of the minimum surface density of solid 
material in the disc can be deduced from the mass and lo- 
cation of the present planets in the solar system: 
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The discontinuity at 2.7 AU stems from the location of the 
ice condensation point (or snow line) that allows the forma- 
tion of icy grains in the outer solar system. Furthermore, 
the total surface density is estimated through the chemical 
composition of the disc, which gives the ratio of gas to solids. 
A fiducial value is 1:0.017 (see Cameron 1973). The surface 
density of the gas component is therefore: 

-3/2 e 
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Since the dust content is rather low, the gaseous component 
is transparent to the visible solar radiation. Thus the gas 
temperature follows from the radiation balance: 
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L is the solar luminosity during the early stages, normalised 
by the present value L & . The three-dimensional density 
structure is given by an isothermal profile 
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with a radially changing scale height h: 
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c s is the sound velocity of an ideal gas with a mean molec- 
ular weight n in units of the hydrogen mass inn- Since the 
density profile is related to a radially varying pressure, the 
gas velocity deviates from the local Keplerian velocity. The 
balance of forces relates the angular velocity fi 9 to the pres- 
sure gradient: 



p dr 



(5) 
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Thus the angular velocity fi 9 of the gas is (see e.g. Adachi 
et al. 1976): 



n g = nv/l-2»fc(r) 

1 dln(p gas cj) / c s V 
% ~ 2 dln(r) l^jsj 



(6) 



It is more appropriate to formulate the rotation of the 
gaseous disc in terms of a velocity lag Av g normalised to 
the local Keplerian velocity Vk- 



Av g = r(Q g - Q) 

~ -VgVK 



(7) 



A typical value of Av g for the minimum mass solar nebula 



at 1 AU is Av a 



-60 m/s. 



This simple model provides a brief description of the 
initial disc. However, the subsequent evolution further mod- 
ifies the structure of the protoplanetary disc. Since embed- 
ded dust grains are coupled to the gas, it is likely that a 
global migration of solid material changes the surface den- 
sity. Moreover, the dust grains are chemically processed, de- 
pending on the local temperature and composition which in- 
troduces additional spatial inhomogeneities. When the grow- 
ing particles pass the critical size of ~ 1 metre, the strong 
onset of radial migration may lead to a final reshaping of 
the distribution of solid material. While these restrictions 
weaken the validity of this approach as the "true" initial 
model, it is still a robust guideline to choose reasonable sur- 
face densities for the solid and the gaseous component after 
the formation of planetesimals. 
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a is the semimajor axis, e is the eccentricity and i is the 
inclination of the orbit. As long as no dominant body is 
structuring the protoplanetary disc, it is justified to assume 
axisymmetry. Hence the argument of the perihelion u, the 
longitude of the ascending node and the eccentric anomaly 
4>e are omitted in the statistical description. 

The deviation of planetesimal orbits from a circle is 
quite small. Thus it is appropriate to expand the above set 
of equations. A planetesimal at a distance ro, in a distance 
z above the midplane and with a velocity (iv, with 
respect to the local circular velocity vk has orbital elements 
(leading order only): 



ro + 2 
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These expressions allow us later a convenient transforma- 
tion between the statistical representation through orbital 
elements and the utilisation of a velocity distribution func- 
tion. 



3 KEPLER ORBITS 

Planetesimals in a protoplanetary disc are subjected to var- 
ious perturbations: Close encounters change their orbits, a 
small but steady gas drag gives rise to a radial drift and 
accretion changes the mass of the planetesimals. While all 
theses processes drive the disc evolution on a timescale of at 
least a few thousand years, each planetesimal moves most of 
the time on an orbit close to an unperturbed Kepler ellipse. 
Though the protoplanetary disc introduces additional per- 
turbations, the central potential dominates for typical disc 
masses around 0.01 Mq. Therefore the classical orbital ele- 
ments still provide a proper framework to study planetesimal 
dynamics. 

The orbital elements of a test particle moving around a 
mass M are: 

_ u 2 _ GM_ 
GM 

OOB(t) = ^ (8) 



4 HILL'S PROBLEM 

When two planetesimals pass close by each other, they ex- 
change energy and angular momentum and separate with 
modified orbital elements. Successive encounters transfer en- 
ergy between planetesimals with different masses, driving an 
evolution of the overall velocity distribution. 

It seems that an encounter is a two-body problem, as 
there are only two planetesimals involved, but the central 
mass has also a major influence turning the problem into 
a three-body encounter The complexity of the problem is 
considerably simplified by reducing it to Hill's problem (Hill 
1878). 

Consider two masses mi and m,2 that orbit a much 
larger mass M c , where both masses are small compared to 
the central mass M c . The mass ratio mi : mi could be arbi- 
trary. This special type of a three body problem is denoted 
as Hill's problem, originally devised to calculate the orbit 
of the moon. It provides a convenient framework to exam- 
ine planetesimal encounters in the potential of a star. The 
equations of motion 4 of the two planetesimals including the 



3 The term three-body encounter does not imply a close passage 
of all involved bodies, but emphasises the strong influence of a 
third one. 

4 The following derivation is quite common to the literature (see 
e.g. Ida 1990; Henon 1986). The later work uses a slightly different 
scaling. 
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central potential and their mutual interaction are: 
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We now introduce the relative vector r and the centre-of- 
mass R: 

ra 2 r 2 + m 1 r 1 



r 2 - ri 
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Furthermore, the equations of motion are transformed to a 
corotating set of coordinates which are scaled by the mutual 
Hill radius rHm and the local Kepler frequency Q 



y = 



rmu 



rmn 
ao((f) — fii) 
rmu 

z' 
rmw 

I mi + m 2 



do 



3M C 



(15) 



where (r',<f>,z') are heliocentric cylindrical coordinates, ao 
is the radius of a properly chosen reference orbit, while the 
mutual Hill radius rmu is an intrinsic length scale of the 
problem : Since the two orbiting planetesimals are small 
compared to the central star, the Hill radius is much smaller 
than the size ao of the reference orbit. Hence it is possible 
to expand the central potential about the reference orbit. 
This yields approximate equations for the relative motion 
and the centre-of-mass: 



x = 2y + 3a; — 3x/r 
y = — 2x — 3y/r 3 
z = —z — 3z/r 3 



X = 2Y + 3X 



(16) 



Y 



-2X 



(17) 



The equations 16 and 17 have some interesting proper- 
ties: Firstly, the centre-of-mass motion separates from the 
interaction of the two bodies. Secondly, the scaled relative 
motion is independent of the masses mi and m 2 , imply- 
ing a fundamental similarity of planetesimal encounters 5 . 
As Eq. 17 is a simple linear differential equation, one read- 
ily obtains the solution of the centre-of-mass motion 



X — b — e cos(t — t) 



Y = -^6t + 7/>- 
Z — i sin(t — ui) 



2e sm(t — t) 



(18) 



which is equivalent to a first-order expansion of a Kepler 
ellipse. u> and r are the longitudes of the ascending node 
and the pericentre, while e and i are the eccentricity and 

5 section 16 makes extensive use of this property. 




Figure 1. Equipotential lines for the effective potential U at 
z = (see Eq. 20). U = U* refers to the largest allowed volume, 
which is enclosed by the Hill sphere and the two Lagrange points 
Li and L 2 . 



inclination scaled by the reduced (i. e. dimensionless) Hill 
radius rmu/ao- The value of b depends on the choice of the 
reference orbit, but it is natural to set b — which implies 
that the centre-of-mass defines the reference orbit. 

While the nonlinear nature of the relative motion (see 
Eq. 16) prevents any general analytical solution, Eq. 18 pro- 
vides at least an asymptotic solution for a large separation of 
the planetesimals, where b can be interpreted as an impact 
parameter. Nevertheless, small b do not necessarily imply 
close encounters, as opposed to the standard definition of 
the impact parameter. However, the expression b — a 2 — ai 
provides a measure of the distance of the two colliding bod- 
ies without invoking the complicated encounter geometry. 
A special solution to Eq. 16 are the unstable equilibrium 
points Li and L 2 at (x,y,z) — (±1,0,0), denoted as La- 
grange points 6 . In addition, an inspection of Eq. 16 reveals 
that the Jacobi energy Ej is conserved: 



■2 , 2 

z + z 



3x 



(19) 



Since the kinetic energy is always a positive quantity, the 
following inequality holds: 



Ej > U = ~ (z 2 



3x 



(20) 



Thus the allowed domain of the particle motion is enclosed 
by the equipotential surfaces of the effective potential U. A 
subset of these equipotential surfaces restricts the allowed 
domain to the vicinity of the origin (see Fig. 1). The largest 
of these surfaces passes through the Lagrange points L\ and 
L 2 . Hence we identify the Hill radius (or Hill sphere) as the 
maximum separation which allows the bound motion of two 
planetesimals 7 . 



6 The additional Lagrange points L3-L5 are missing due to the 
Hill approximation. 

7 The same argument applies to the tidal boundary in cluster 
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Beside the numerical solution of the equations of mo- 
tion, it is useful to define osculating (or instantaneous) or- 
bital elements 



b = 4x + 2y 

2 2,-2 
— Z + Z 

2 = x 2 + (3x + 2yf 



(21) 



t[> = -2x + y + -bt 

uj = t — arctan(0, z) 

t = t — arctan(£, 3a; + 2y) 



1 / 2 . -2-, 3,2 

-(e + i ) - -b 



(22) 



(23) 



which provide a convenient description of the initial relative 
orbit and the modified orbit after the encounter. 



5 PROTOPLANET GROWTH 

Our work follows the evolution of a planetesimal disc into 
few protoplanets, including the full set of interaction pro- 
cesses. Hence we summarise the main aspects of protoplanet 
growth first to provide a robust framework. 

Although the sizes of the planetesimal cover a wide 
range, they virtually form two main groups: The smaller 
field planetesimals and the embedded protoplanets (or their 
precursors). This two-group approximation (e.g. Wetherill 
1989; Ida 1993) allows one to have a clearer insight into the 
growth process. 

During the initial phase all planetesimals share the same 
velocity dispersion independently of their mass. The initial 
random velocities are low enough for an efficient gravita- 
tional focusing. Hence, the growth rate of a protoplanet with 
mass M radius R can be estimated as (e.g. Ida et al. 1993): 



M ■ 



E B 2 



1 + 



2GM \ 



(24) 



Eq. 24 is the two-body accretion rate, which should be mod- 
ified due to the gravity of the central star. Nevertheless this 
approximation gives an appropriate description to discuss 
the basic properties of the growth mode. The scale height 
H (Eq. 215) and the relative velocity u re i are related to the 
mean eccentricity e m = sj (e 2 ) of the field planetesimals: 



H « w r ei/n «rd ~ e m aQ, 



(25) 



Thus the accretion rate in the limit of strong gravitational 
focusing (2GM/R > v 2 cl ) is: 

EGM 



M « 2ttR 



a 2 fie 2 



M 4/3 



(26) 
(27) 



If the protoplanets are massive enough, they start to control 
the velocity dispersion of the planetesimals in their vicinity. 



dynamics or the Roche lobe in stellar dynamics, which are equiv- 
alent to the Hill sphere. 



The width Aa of this sphere of influence, the heating zone, 
is related to the Hill radius Rum of the protoplanet (Ida 
1993): 



(28) 



Aa = AaR mll = 4R H iin/ ^{e 2 m + i 2 m ) + 12 



h = 



M 
3M~ C 



e m and i are eccentricity and inclination of the field plan- 
etesimals, scaled by the reduced Hill radius h of the pro- 
toplanet. M c is the mass of the central star. The condition 
that the protoplanet controls the velocity dispersion of the 
field planetesimals reads (Ida et al. 1993): 



2M 2 



> Em 



(29) 



27raAa 

This condition is equivalent to a lower limit of the proto 
planetary mass: 

3 / 5 / V„2 \ 3/5 / n 3/5 



M 



> 



ivAaY 



Sa 2 

M c 



m 



(30) 



M/m depends on several parameters, but reasonable values 
yield M/m w 50-100. The velocity dispersion in the heated 
region is roughly 



v w Rmnti 



(31) 



which gives an interesting relation to the condition that 
leads to gap formation. A protoplanet can open a gap in 
the planetesimal component if it is larger than a critical 
mass M gap (Rafikov 2001) 



E 
Al, 



M c 



E 
M, 



V3 
1/3 



(rn_ V 
lc \M C ) 

lc \ M c ) \ v J 



if v < Qr Hiu 
if v » firHiii 



(32) 



where rmn is the Hill radius of the field planetesimals. If the 
velocity dispersion v is controlled by the protoplanet, Eq. 32 
together with Eq. 31 demonstrate that the condition for gap 
formation is equivalent to Eq. 30, i. e. the efficient heating of 
the field planetesimals implies gap formation and vice versa. 
The higher velocity dispersion of the field planetesimals (see 
Eq. 31) reduces the growth rate given by equation 26 to 



.RRmn 

52 



ocM 2/3 



(33) 
(34) 



Different mass accretion rates imply different growth mode. 
If two protoplanets have different masses Mi and M 2 , their 
mass ratio evolves as: 



_d Ma 
dt MT 



Ma 
Mi 



Ma 
M 2 



Mi 
Mi 



(35) 



When the growth timescale M/M decreases with mass, a 
small mass difference increases with time. This is the case 
for Eq. 26, which gives rise to runaway accretion. As soon as 
the protoplanets control the velocity dispersion of the field 
planetesimals, the growth timescale increases with mass and 
therefore the protoplanet masses become more similar as 
they grow oligarchically. 

The field planetesimals also damp the excitation due 
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E[g/cm 2 ] 


M iso /M Q 


M iso /M ffi 


2 


3.91 x 10~ 8 


0.013 


10 


4.33 x 10~ 7 


0.144 


100 


1.37 x 10~ 5 


4.548 



timescale 



Table 1. Isolation mass for different surface densities at r = 1 
AU and M c = 1% 



to protoplanet-protoplanet interactions and keep them on 
nearly circular orbits. The balance between these scatterings 
and the dynamical friction due to smaller bodies establishes 
a roughly constant orbital separation b (Kokubo 1997): 



b = R f 



5/ 7e 2 m M 
2-rvEaR 



b = b/Rmn 



(36) 
(37) 



R is the radius of the protoplanet, M is its mass and e 
is the reduced eccentricity of the field planetesimals. The 
stabilised spacing prevents collisions between protoplanets, 
but it also restricts the feeding zone - the area from which a 
protoplanet accretes. If all matter in the feeding zone is ac- 
creted by the protoplanet, it reaches its final isolation mass 
(Kokubo & Ida 2000): 



M iso = 2-nbaE 



(38) 



Inserting Eq. 36 yields the isolation mass in units of the 
mass of the host star M c : 

M iso /M c = (112^) 3/8 ( \Y" f?V /8 (e-) 3/8 



a 2 E\ 3/2 /aV 1/8 



M c J \M C 
M c Mc 



19.67 x (e 2 m ) 



2 \3/8 



(39) 



There is a weak dependence on the density p of the pro- 
toplanet, but the most important parameter is the surface 
density E. If we take the minimum mass solar nebula as an 
example, the radial dependence of the surface density im- 
plies an isolation mass that grows with increasing distance 
to the host start. Hence protoplanets beyond some critical 
radial distance are massive enough (larger than « 15 
Bodenheimer 1986) to initiate gas accretion from the proto- 
planetary disc. 

As the protoplanets approach the isolation mass, inter- 
actions with the gaseous disc and neighbouring protoplanets 
become increasingly important. We estimate the onset of or- 
bit crossing by a comparison of the perturbation timescale 
Tpcrt of protoplanet-protoplanet interactions with the damp- 
ing timescale Tdamp due to planetesimal-protoplanet scatter- 
ings. Since the protoplanets are well separated (6 w 5 . . . 10), 
it is possible to apply perturbation theory (Petit 1986,see 
e.g.): 



^~pcrt 



b b 

7hQ, 



(40) 



''"damp 



n3/2 



2ttG 2 ln(A)(M + m)n m 



(41) 



where T r and T z are the radial and vertical velocity disper- 
sion of the field planetesimals. Hence the criterion for the 
onset of orbital crossing is: 



Tpert *C Tdamp 



(42) 



As the protoplanets control the velocity dispersion of the 
field planetesimals (see Eq. 31), this condition reduces to: 



E M > E m ln(A) 
> E m x / 



72 
7tt 



(43) 



We anticipate section 16 (see Eq. 221) to derive the damping 



Thus orbital crossing sets in when the mean surface den- 
sity Em of the protoplanets exceeds some fraction / of the 
field planetesimal density E m . While the factor / depends 
strongly on the separation b of the protoplanets, a fiducial 
value is / « 1, in agreement with the estimates of Goldreich 
et al. (2004). 

The onset of migration and the resonant interaction 
of protoplanets with the disc and other protoplanets ter- 
minates the local nature of the protoplanet accretion pro- 
cess and requires a global evolution of the planetary system. 
While the final stage deserves a careful analysis, further re- 
search is beyond the scope of this work. 



6 DIRECT-SUMMATION TECHNIQUES 

FROM THE STANDPOINT OF PLANETARY 
DYNAMICS: FIRST STEPS 

The protoplanet formation is essentially an A-body prob- 
lem. Although we seek for a more elaborated solution to this 
problem which benefits from statistical methods, the pure 
A-body approach is a logical starting point. Direct calcu- 
lations with a few thousand bodies have provided us with 
a valuable insight into the growth mode (Ida 1992; Kokubo 
1996, see e.g.), but they are also powerful guidelines that help 
developing other techniques. Statistical calculations rely on 
a number of approximations and "exact" A-body calcula- 
tions provide the necessary, unbiased validation of the de- 
rived formula. 

The choice of the integrator is a key element in the 
numerical solution of the equations of motion. Our require- 
ments are the stable long-term integration of a few ten thou- 
sand planetesimals with the capability of treating close en- 
counters, collisions and the perspective to evolve it into an 
improved hybrid code. Approximative methods like the Fast 
Multipole Method or Tree codes have a scaling of the com- 
putational time close to N, but the accuracy in this regime 
is too poor to guarantee the stable integration of Keple- 
rian orbits (compare the discussion in Hernquist et al. 1993; 
Spurzem 1999). 

The class of exact methods (all scale asymptotically 
with N 2 ) roughly divides in two parts: 

(i) Symplectic methods (see e.g. Wisdom 1991 or the 
Symba code, Duncan et al. 1998) rely on a careful expan- 
sion of the Hamiltonian which guarantees that the numerical 
integration follows a perturbed Hamiltonian. While there is 
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still an integration error, all properties of a Hamiltonian sys- 
tem like conservation of phase-space volume are conserved 
by the numerical integration. The drawback of these very el- 
egant methods is that the symplecticity is immediately bro- 
ken by adaptive time steps, collisions or complicated exter- 
nal forces if no special precaution is taken. Even the numeri- 
cal truncation error breaks the symplecticity to some extend 
(Skeel 1999). Moore & Quillen (2010) recently developed a 
parallel integrator for graphics processing unit (GPU) that 
uses symplectic and Hermite algorithms according to the 
resolution needed. The algorithm is similar to Symba, but 
less accurate. 

(ii) The second group represents the "classical" methods 
that are based on Taylor expansions of the solution. They 
come in different flavours like implicit methods, predictor- 
corrector integrators or iterated schemes. Time symmetric 
methods stand out among these different approaches, as 
they show no secular drift in the energy error. These in- 
tegrators are so well-behaved that one may even call them 
"nearly symplectic" . 

Taking all the requirements into account we have chosen 
Nbody6++ 8 , an integrator which is the most recent descen- 
dant from the famous TV-body code family from S. Aarseths' 
"factory". This version was parallelised by Spurzem (1999), 
which opened the use of current supercomputers. 

This parallel version, named Nbody6+ + , offers many 
versatile features that were included over the past years and 
more and more refined as time passed by. While all these 
elaborated tools deserve attention, we restrict ourselves for 
brevity to the components which contribute to the planetes- 
imal problem. The main components of the code are: 

(i) Individual time steps and a block time step scheme. 

(ii) Ahmad-Cohen neighbour scheme. 

(iii) Hermite scheme. 

(iv) KS-Regularisation for close two-body encounters. 

We will explain each of these features and highlight their 
advantages for the main goal of this work. This is impor- 
tant to understand the immediate first level of modifica- 
tions required to adapt the numerical scheme to planetary 
dynamics. We present later, in section 8, the next level of 
complexity in the modifications to be carried out to mold 
the numerical tool to our purposes. 

6.1 Individual Time Steps 

The choice of the time step controls the accuracy as well 
as the efficiency of any given integrator. Too small time 
steps slow down the integration without necessity, whereas 
too large values increase the error. An efficient solution to 
this dilemma is an adaptive time step that is adjusted af- 
ter each integration step according to a specified accuracy 
limit. While the idea is quite clear, there is no unique receipt 
how to choose the proper time step. A common approach for 
iV-body systems is to use the parameters at hand (like par- 
ticle velocity, force, etc.) to derive a timescale of the particle 
motion. This procedure leaves enough space for a wealth of 

8 Aarseth (1999) gives a nice review on the remarkable history 
of the NBODY-codes, more details are given in Aarseth (2003). 



different time step criteria. Nbody6++ uses the standard 
Aarseth expression (Aarseth 1985) 

I FF( 2 ) + (F(!)) 2 

At = ^ FWF , 3) + {Fi2))2 

which makes use of the force and the time derivatives up to 
third order. 

The time step choice is not unique in multi-particle sys- 
tems. One solution is to take the minimum of all these values 
as a shared (or global) time step, but this is not recom- 
mended unless all individual steps are of the same size. 

The second option is to evolve each particle track with 
its own, individual time step. This method abandons the 
convenient synchronisation provided by a global step, there- 
fore each integration of a particle demands a synchronisation 
of all particles through predictions. Since the prediction of 
all particles is O(N), it is counter-balanced by the saving 
of force computations. Nevertheless the overhead is still sig- 
nificant, so a further optimisation might be desired. The 
basic idea of the block time step method is to force parti- 
cles in groups that are integrated together, which reduces 
the number of necessary predictions by a factor comparable 
to the mean group size. These groups are enforced through 
two constraints. The first condition is a discretisation of the 
steps in powers of two: 

At = 2~ k k^O (45) 

This condition increases the chances that two different par- 
ticles share the same timelevel, but it also reduces roundoff 
errors since the time steps are now exactly representable 
numbers. The second condition locks the "phases" of parti- 
cles with the same time step 

Ti = mod AU (46) 

i. e. the particle time T» is an integer multiple of the actual 
time step AU so that all particles with the same step share 
the same block. Note that a step can not be increased at 
any time T;, but only when the second condition (Eq. 46) 
for the larger step is met. 

6.2 Ahmad— Cohen Neighbour Scheme 

The first TV-body codes calculated always the full force (i. e. 
summation over all particles) to integrate a particle. But 
not all particles contribute with the same weight to the to- 
tal force. Distant particles provide a smooth, slowly vary- 
ing regular force, whereas the neighbouring particles form a 
rapidly changing environment which gives rise to an irregu- 
lar force. The Ahmad-Cohen neighbour scheme (Ahmad & 
Cohen 1973) takes advantage of this spatial hierarchy by di- 
viding the surrounding particles in the already mentioned 
two groups according to the neighbour sphere radius R s - 
Both partial forces fluctuate on different timescales, which 
are calculated according to Eq. 44. The key to the efficiency 
of the method is the inequality 

Atirr < At rcg (47) 

Regular forces are extrapolated between two full force cal- 
culations 

F rog = F£> + At Fill + \ (A*) 2 FL 2 g + \ (At) 3 F< c 3 g (48) 
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Figure 2. All particles inside the neighbour sphere R a are se- 
lected as neighbours (filled circles). As the neighbour list is fixed 
during regular steps (marked by black lines), a second shell in- 
cludes possible intruders. 



with a third order accurate expression, whereas the plain 
NbodyGH — h uses only linear extrapolation (see the next sec- 
tions for a detailed discussion). 



6.3 Hermite Scheme 

The Hermite scheme is a fourth-order accurate integrator 
that was applied first by Makino (1992) to the integra- 
tion of a planetesimal system. This scheme is used to in- 
tegrate single particles and CM-bodies 9 in the main part of 
Nbody6H — K It is accomplished by a predictor and a correc- 
tor step. The prediction is second order accurate: 



x p = x + v At - 
v p = v + a At - 



ia (At) 2 + ^(At) 3 
Ja (At) 2 



(49) 



Now the acceleration is evaluated at the predicted position 
to derive the second and third order time derivatives of the 
force: 



,(3) 



-4a - 2a p 

At 
6ao + 6a p 



— 6ao + 6a p 

At 2 

12ao — 12a p 



(50) 

At 2 + A t3 "~* (51) 
Using the additional derivatives one can improve the predic- 
tion: 

x c = x p + ^a< 2) (At) 4 + ^^(At) 5 + 0(At 6 ) 



20 



(At) 3 + C>(At 6 



20 



7 . 1 . 

— an H a, 

60 30 1 



(52) 



v, + ia^(At) 3 + iaf (At) 4 



v p + - (-a + a p ) At + 



12 



an 



0(Ar 
l . 

-12 aj 



C(At J 



(At) 2 + 
(53) 



9 CM denotes centre-of-mass, see the section on KS- 
Regularisation for more details. 



The corrected positions are fourth order accurate. While the 
Hermite scheme is robust and stable, even in combination 
with the neighbour scheme, it is not accurate enough to 
integrate planetesimal orbits efficiently. 

6.4 Hermite Iteration 

The Hermite scheme bears the potential for a powerful ex- 
tension, since it is a predictor-corrector scheme. An essential 
part of this scheme is the calculation of the new forces with 
the predicted positions, but it should improve the accuracy 
if they are recalculated using the corrected positions. The 
new forces are readily used to improve the corrected values, 
which closes the scheme to an iteration loop - the Hermite 
iteration (Kokubo et al. 1998). 

There are only few modifications necessary to obtain 
the iterated scheme. The particle prediction remains second- 
order accurate: 

1 9 1 t 

Xp = x + v At + -an(At) 2 + -a (Atr 



vo + a Af + -an (At) 2 



(54) 



The force and its first time derivative are calculated to derive 
higher derivatives according to Eq. 50 and 51: 

a p = f(x p , v p ) (55) 
,v„) (56) 



a p — f(x p , v,, j 



The new corrector omits the highest order term in the posi- 
tion, making the velocity and the position to the same order 
accurate: 

x c = Xp + ^a 2) (At) 4 + O(At 5 ) 



= Xp + -(a p - a )(At) + 
(At) 3 + C>(At 5 ) 



Tan 



12 



1/ / 5. 1 . 

2 (-an +a p )At+ I -^a - — a p 



(At) 2 + 0(At 5 ) 



(57) 



(58) 



It seems unreasonable to drop one order in the position, but 
a reformulation of the predictor-corrector step reveals that 
this slight change yields a time symmetric scheme: 



1 1 2 

v c = v + -(a, + an) At - T^(a p - an) (At) 

1 1 2 

x c = xn + -(v c + v )At - ^(ap - an) (At) 



(59) 



The iteration is achieved by returning to Eq. 55-56 with the 
predicted positions replaced by the more accurate values 
x c ,v c . Although the integration does not need the second 
and third time derivative of the forces explicitly, it is useful 
to provide them at the end of the iteration for the calculation 
of the new time step: 

2an + 4a p 6an — 6a p 



a (2) (t + At) 
a (3) (t + At) 



At 
6ao + 6a p 



At 2 
12an — 12a p 



(60) 

At 2 ' •~"a,^ ( 61 ) 
Numerical tests show that convergence is reached after one 
or two iterations, making this approach very efficient. It 
needs less force evaluations than the Hermite scheme for 



Hybrid methods: A new composite algorithm 11 



the same accuracy. Moreover, its time symmetry suppresses 
a secular drift of the energy error. 



6.5 Extended Hermite Scheme 

The Hermite scheme is an integral part of Nbody6+ + 
and proved its value in many applications. It would have 
been natural to improve the performance with an additional 
iteration, but our first tentative implementations showed 
rather negative results: The iterated scheme was more un- 
stable, slower and even less accurate than the plain Hermite 
scheme. An inspection of the code structure revealed that 
the Ahmad-Cohen neighbour scheme is the cause of this. 

Each particle integration is composed of two parts - 
frequent neighbour force evaluations and less frequent total 
force evaluations including derivative corrections. Every reg- 
ular correction leads to an additional change in the position 
of a particle, which introduces a spurious discontinuity in 
the neighbour force and its derivatives. The Hermite itera- 
tion reacts to this artificial jump in two ways: It increases the 
regular correction, and - what is more important - it am- 
plifies any spurious error during the iteration which leads to 
an extreme unstable behaviour. 

Since the Hermite iteration is a key element in the ef- 
ficient integration of planetesimal orbits, we sought for a 
modification of the Hermite scheme that circumvents the de- 
picted instability. The problem gives already an indication of 
a possible solution. A scheme with much smaller corrections 
would not suffer from the feedback of spurious errors. 

Nbody6++ stores already the second and third time 
derivative of the forces for the time step calculation. A man- 
ifest application of these derivatives at hand is the improve- 
ment of the predictions to fourth order: 



1 o 1 q 

x p = x + v Ai + 2 a o(At) + -ao(At)' 3 



, = v + a At+ ^ao(At) 2 
+ iaf(A^ + la ( f(At) 4 



(62) 



(63) 



The prediction to fourth order was used in the iterative 
schemes of Kokubo et al. (1998); Mikkola & Aarseth (1998). 
Again, the new forces a p and a p are calculated to improve 
but with a modified corrector: 



Finally, the derivatives are updated: 

& (t + At) = a p 
ao(t + At) = a p 
ao (2) (t + At)=a( 2 >+Atai 3 > 



a (3 \t + At) = a 



(3) 



(70) 
(71) 
(72) 
(73) 



The new scheme has an appealing property, which is related 
to the usage of the higher force derivatives. As the predictor 
is fourth-order accurate, it is equivalent to one full Hermite 
step. Since the corrector uses new forces to improve the two 
highest orders, it is equivalent to a first iteration step. Thus 
we obtained a one-fold iterated Hermite scheme at no ex- 
tra cost. This extended Hermite scheme reduces the energy 
error by three orders of magnitude, compared to the plain 
Nbody6++ with the same number of force evaluations. 



6.6 KS— Regular isat ion 

Two bodies undergoing a close encounter are integrated in 
a special set of regular coordinates that separates the rela- 
tive motion from the motion of the centre-of-mass. A close 
encounter poses a numerical problem due to the singular- 
ity of the gravitational forces at zero separation. While the 
growing force amplifies roundoff errors as the two bodies 
approach each other closely, the collision is only an appar- 
ent singularity since the analytic solution stays well-defined. 
This opens the possibility of a proper coordinate transfor- 
mation which removes the singularity from the equations of 
motion. The Kustaanheimo-Stiefel regularisation takes ad- 
vantage of a four-dimensional set of variables to transform 
the Kepler problem into a harmonic oscillator (Kustaan- 
heimo & Stiefel 1965). Perturbations are readily included 
in the new set of equations of motion. 

The centre-of-mass is added as a pseudo-particle, the 
CM-body, which is integrated as a normal particle plus a 
perturbation force due to the deviation from a point mass. 
See Mikkola (1997) or Mikkola & Aarseth (1998) for more 
details. 



x p and v p 



ap = f(xp,v p ) (64) 7 ADDITIONAL FORCES FOR 

a p = f ( X p Vp) (65) PLANETESIMAL DISC DYNAMICS 

(2) . . — 4ao — 2a p — 6ao + 6a p . . Nbody6H — h only includes the gravitational interaction of all 
a ™ ' ' At ^ At' 2 particles, therefore additional forces have to be added "per 

(3) , . 6ao + 6a p 12ao — 12a p , . hand". A planetesimal disc requires two new forces: The 
n ^ ' At 2 At 3 presence of a central star introduces an additional central 

1 . ( 2 ) (2)\/ A \4 potential, while the gaseous component of the protoplane- 

Xc Xp 24 " a ° taxy disc is the source of a friction force. It is important 

1 , (3) ( 3 )wArt 5 ffisl that the new forces are properly included in the neighbour 

120 a ° scheme to assure that regular steps remain larger than ir- 

v _ v 1 / a (2) a ( 2 h(At) 3 -|- regular steps. Since a dissipative force breaks the energy 

p 6 conservation, one has to integrate the energy loss as well to 

1 / (3) (3)wa.\4 maintain a valid energy error control. In the next subsections 

— {a. K n '~a yl )(At) (69) ° J 

24 we describe how we have done this. 
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7.1 Central Potential 

A star is much heavier than a planetesimal. Thus, the central 
star is introduced as a spatially fixed Keplerian potential: 



<E> C 



GM C 



GmM c x 



(74) 



Since the orbital motion of the planetesimals sets the domi- 
nant (and largest) dynamical timescale in the system, we in- 
cluded the central force as a component of the regular force. 
Moreover, the central potential also introduces a strong syn- 
chronisation, since planetesimals in a narrow ring share vir- 
tually the same regular block time step. 



7.2 Drag Force 

As the whole planetesimal system is embedded in a dilute 
gaseous disc, each planetesimal is subjected to a small, but 
noticeable drag force. The drag regime 10 depends on the gas 
density and the size of the planetesimals. Kilometer-sized 
planetesimals are subjected to the deceleration 



dv 

dt 2m 
C D = 0.5 



'"pgas-RV - V 9 |(V - V 9 ) 



(75) 
(76) 



which is inversely proportional to the radius R(m) of the 
planetesimal in this drag regime. v g is the rotational velocity 
of the gaseous disc, which rotates slower than the planetes- 
imal system as it is partially pressure supported. The drag 
force leads to an orbital decay a of the semimajor axis of a 
planetesimals: 



3 gj _ Pgas {(AV) 2 } 

4 PBody R(m)Q 



(77) 



Thus smaller particles migrate faster, with a maximum at 
R « 1 m. Even smaller bodies couple to the gas, which 
reduces the effective drag force. The dissipation rate and its 
time derivative are: 



Wdrag — Fdrag ' V 

Wdrag = Fdrag ' V + Fdrag " F to t 



(78) 



We integrate the dissipation rate Wdmg to maintain a valid 
energy error: 



AE : 



/•t i 



Wdrag dt 



At = t 2 - ti 

AE = i (VKdrag.l + Wdrag, 2 ) At 



+ 



hi** 



1 - Wdrag,2 At + Q(At° 



(79) 
(80) 

(81) 



The expression is fourth order accurate in accordance with 
the order of the extended Hermite scheme. 



7.3 Accurate integration of close encounters: 

Tidal perturbations of KS— Pairs and impact 
of the gaseous disc 

Both new forces also demand a modification of the regu- 
larisation treatment. They perturb the relative motion of a 
KS-pair and modify the orbit of the centre-of-mass. While 
the modification of the equations of motion is rather clear, 
the neighbour scheme requires some additional work. 
Let ri , r 2 be the positions of the two regularised particles. 
The equations of motion read (G = 1) 

r = r2 — ri 
ri 



ri = -Afc-g + rn 2 ^ + F drag ,i 
r 1 r 

V2 r 

r 2 = -M C ^T - mi — + Fdrag, 2 



(82) 



where the perturbations by other particles have been omit- 
ted for clarity. Centre-of-mass motion and the orbital motion 
are separated: 



r = - M c ^| + M c ^ + Fdrag,2 - Fdrag.l 

r r 2 r l 

f? _ Af miri M 7712 1-2 ->- mi V 
M r{ M r| M 



mi 
~M 



Fdr 



r = r2 — ri M — mi + m.2 



(83) 

(84) 
(85) 
(86) 



Two new contributions show up due to the external 
forces: The KS-pair is tidally perturbed by the central star 
and influenced by the gaseous disc. While the aerodynamic 
properties of a single particle are well understood, two bod- 
ies revolving about each other may induce complex gas flows 
in their vicinity, which could invalidate the linear combina- 
tion of the drag forces on each component. Therefore we drop 
the drag force term to avoid spurious dissipation. Since the 
dynamic environment allows virtually no stable binaries 11 
in a planetesimal disc, the influence of the drag force on the 
encounter dynamics is negligible. 

We further decompose the additional acceleration of the 
centre-of-mass motion, since the neighbour scheme benefits 
from a clear separation of the timescales. Therefore, the tidal 
perturbation is split into a smooth mean force and a pertur- 
bation force: 



= -M c 



R. 

R 3 

R mi ri 

: i?3 c Mr\ 

= + O(r 2 ) (88) 
The mean forces varies on the orbital timescale and is hence 



JTl2 V2 



F pcrt =M c ^-M c ^-M c ^ (87) 



10 The main drag regimes are Stokes (laminar flow), Epstein 
(mean free molecular path larger than object size) and Newton's 
drag law (turbulent flow). Wcidcnschilling (1977) provides a nice 
review on the different drag regimes. 



Tidal capturing of moons starts in the late stages of planet 
formation, but is limited to the planets or their precursors. How- 
ever, the quiescent conditions in an early Kuiper belt lead to a 
more prominent role of binaries. See the summary of Astakhov 
ct al. (2005). 
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Time Step Ratio 



1000 



i>= 100 




r/i'H 



Figure 3. Time step ratio for 7V n b = 100. Curves are plotted for 
different values of a v /(rumQ). The dotted line is approximation 
Eq. 95. 



included as a regular force component, while the perturba- 
tion is treated as an irregular force as it changes with the 
internal orbital period of the pair. 



however, that the relation of regular and irregular costs is 
more complicated with GPU technology) 



At- 



FFW + (F(!)) 2 



(2)12 



(89) 



where F^ are the force and its time derivatives. It is applied 
to the calculation of the regular step using the regular force 
and accordingly to the irregular step based on the irregular 
force. The regular time step of a particle orbiting the central 
mass Me at a distance ro is: 



At r = 



tT)r 



Q = 



GM C 



(90) 



For simplicity, we introduce the scaled timestep ratio jt = 
JtyVbr/Vtee- The free parameters of the problem are the 
mean particle distance f, the velocity dispersion a v (addi- 
tional to the Keplerian shear), the particle mass m» and the 
neighbour number iV n b. We employ Hill's approximation for 
the central potential and obtain: 



8 FURTHER TUNING THE DIAL: 

OPTIMISING TV BODY TO THE DISC 

An astrophysical simulation is a tool to analyse problems 
and predict dynamical systems which are not accessible 
to experiments. The design of a new simulation tool does 
not only require the careful implementation of the invoked 
physics, but also an analysis of the code performance to 
make best use of the available hardware. 

We want to apply Nbody6+ + for the first time to 
planet formation, a subject that is quite different to stellar 
clusters. The central star forces the planetesimals on regu- 
lar orbits which need higher accuracy than the motion of 
stars in a cluster. In addition, the orbital motion also in- 
troduces a strong synchronisation among the planetesimals, 
thus allowing a more efficient integration. 

We examine the differences due to the integration of a 
disc system in the following sections. In particular, we will 
address in detail the role of the geometry of the problem 
and the neighbour scheme, the prediction of the number of 
neighbouring particles, the communication, the block size 
distribution and the optimal neighbour particle number for 
the direct-summation of the massive particles in the proto- 
planetary system 

8.1 Disc Geometry and Neighbour Scheme 

The introduction of the neighbour scheme by Ahmad & Co- 
hen (1973) has provided us with a technique to save a con- 
siderable amount of computational time in star cluster simu- 
lations. Since the average ratio of the regular to the irregular 
time step jt is of the order of 10, the integration is speeded 
up by the same factor. One may expect a similar speedup 
for planetesimal systems, but in this case the time step ratio 
is roughly three. The time step is calculated with the stan- 
dard Aarseth time step criterion (it should be mentioned, 



AU » ^/(«, rmu, f , <Tv, N nb ) (91) 

— 3 M ™ 

f is a yet unknown function. Dimensional analysis leads to 
M^^f(^-',N nh ) (93) 

il \ r Hill" T'Hill / 

Vrmiifi rum J 

which shows that the time step ratio is essentially con- 
trolled by the interparticle distance and the velocity disper- 
sion. We generated different random realisations of planetes- 
imals discs with different densities and velocity dispersions 
to cover the range of possible values. The neighbour number 
is fixed to 7V n b = 100 to reduce the noise due to small num- 
ber statistics, but 7 t converges to a value independently of 
the neighbour number already for A n b > 10. Fig. 3 shows 
the numerical calculation of the time step ratio for various 
values of f and a v . A good approximation to the calculated 
values of jt is: 



1.79x^1 + 1.03^+0.94^ (95) 

Planetesimal discs have usually a small velocity dispersion 
(compared to the orbital velocity) and a low density in terms 
of the Hill radius, which leaves a major influence to the Ke- 
plerian shear. Since the shear motion is directly linked to 
the local Keplerian frequency, this synchronisation reduces 
"/t to values smaller than ten. The numerical calculations 
show larger time step ratios with increasing velocity disper- 
sion and for high densities 12 , but planetesimal discs are far 
from these extreme parameter values. 
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(1) V-Exponent=0 

(2) M-Exponent=l/6 



— I— 
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Figure 4. Regular steps per particle and per 1 TV-body time 
in the inner core (Ri < 0.5) of a 5000 particle plummer model. 
Plotted are (1) different mass exponents with velocity exponent 
and (2) different velocity exponents with mass exponent 1/6. 



8.2 Optimal Neighbour Criterion 

The standard neighbour criterion uses the geometrical dis- 
tance: Particles are neighbours if their distance to the refer- 
ence particle is smaller than a limit R s . This criterion is sim- 
ple and probably the best choice for an equal mass system. 
However, a multi-mass system may require a different crite- 
rion, since a massive particle outside the neighbour sphere 
could have a stronger influence than lighter particles inside 
the neighbour sphere. Also, relative effects are smaller at 
large distances. A more appropriate selection should rely on 
some "perturbation strength" of a particle. 

It turned out that a better criterion is the magnitude 
of the fourth time derivative of the pairwise force F^' , 
i. e. those particles are selected as neighbours which pro- 
duce the largest integration error in accordance with the 
Hermite scheme. is a complicated expression (compare 
Appendix A), but the leading term can be estimated via 
dimensional analysis: 



(96) 



We use this expression to define a new apparent distance 
between the integrated particle % and a neighbour j: 

1/6 / \ 2/3 

(97) 

v a is an arbitrary scaling velocity to obtain a distance with 
dimension length. This new distance definition moves mas- 
sive or fast neighbours to an apparently smaller distance, 
thus enforcing that these particles are preferentially included 
in the neighbour list. In addition, the modified distance is 
readily included in the conventional neighbour scheme. We 
tested different mass and velocity exponents to verify that 
Eq. 97 is the optimal choice. Figure 4 shows that these ex- 
ponents are indeed the optimal choice for a Plummer model 
(Plummer 1911) with mass spectrum. The new scheme saves 
25% of the force evaluations in the core, but the impact on 

12 P/fHUl < I corresponds to an unstable self-gravitating disc. 



a planetesimals system is smaller, as it is the case for the 
neighbour scheme. While a velocity dependent distance re- 
duces the number of necessary full force evaluations, it in- 
troduces a distance changing with time which destabilises 
the integration. The result is a much larger energy error 
compared to the achieved speedup. Therefore we only rec- 
ommend the mass modification of the apparent distance. 



8.3 Neighbour Changes 

The rate at which the neighbours of a given particle change 
has a noticeable influence on the accuracy of the code. Dur- 
ing the course of an integration the second and third time 
derivative of the regular and irregular force are calculated 
from an interpolation formula (see Eq. 50 and Eq. 51). 
Whenever a particle leaves (or enters) the neighbour sphere, 
these derivatives are corrected by analytic expressions 13 . 
Hence many neighbour changes lead to a pronounced spuri- 
ous difference. 

We estimated the rate at which particles cross the 
neighbour sphere boundary to quantify this effect. Neigh- 
bour changes are due to the Keplerian shear and the super- 
imposed random velocities of the particles. The two effects 
lead to: 



N +/- 
N +/- 



= At 



Nnb 
-L orb 



3 . , T &v 
-At r N nb —^- 

2 -Knb 



At r 



N n h 3iva v 

Torb RnbQ 



Shear 



Dispersion 



(98) 



7?nb and N n b are the neighbour sphere radius and the 
neighbour number, respectively. In practice, the neighbour 
changes due to the shear account for up to 80% of the 
total neighbour changes. The standard regular time step 
At r = 2~ 5 and 50 neighbours yield a change of one par- 
ticle per regular step, which is fairly safe. 



8.4 Neighbour Prediction 

Each integration step is preceded by the prediction of all 
neighbours of the particles that are due. A regular step re- 
quires the full prediction of all particles, so there is no possi- 
bility to save computing time. In contrast, an irregular step 
calculates only neighbour forces, which requires the predic- 
tion of less particles. Thus the prediction of all particles 
to prepare an irregular step is a simple, but, depending on 
the block size, computational costly solution. It seems to be 
more efficient to predict only the required particles, but ran- 
dom access to the particle data and the complete check of 
all neighbour list entries introduces an additional overhead. 
Therefore large block sizes should favour the first approach, 
whereas the second approach is more suitable for small block 
sizes. 

Both regimes are separated by a critical block size AT;* r . 
If Ni IT particles with (N^) neighbours are due, then only 
JV merge particles need to be predicted: 



13 Appendix A gives a complete set of the force derivatives up to 
third order. 
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Table 2. Ring Communication. 


Communication partners are 


Irregular 
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0.35 
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877 


8.8 


fixed, while the exchanged data 


varies, rij 




1 cycles are needed. 


Regular 


20 


0.22 


368 


1668 


75 



Table 3. Hierarchical Communication. Communication partners 
change after every cycle. The exchanged data amount doubles 
with every new cycle, hence only ln2(n p ) cycles are needed. 
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a ~ iVtot ^1 - exp ^~ 

iV irr (iVnb) 



N irl {N n 

JVtot 



(99) 



The size -/V mcrge of the merged neighbour lists is smaller 
than the total number of neighbour list entries, since some 
particles are by chance members of more than one neighbour 
list. Performance measurements show that the prediction of 
the merged neighbour lists is 10 % more costly (per particle) 
than the full prediction, mainly due to additional sorting and 
a random memory access. Thus N* IT satisfies: 



1.1 X AT 



Inserting Eq. 99 yields the critical block size: 

Mot 



N* 

irr 



2.4 



<AU> 



(100) 



(101) 



The prediction mode is chosen according to the actual block 
size. 



8.5 Communication Scheme 

Nbody6++ is parallelised using a copy algorithm. A com- 
plete copy of the particle data is located on each node, so 
the integration step of one particle does not need any com- 
munication. Therefore a block of Nbi particles is divided 
in n v parts (n p is the processor number), which are inte- 
grated by different processors in parallel. The integration 
step is completed by an all-to-all communication of the 
different subblocks to synchronise the particle data on all 
nodes. Hence the amount of communicated data is propor- 
tional to JVbi x n p . A communication in a ring-like fashion 
(see table 2) needs n p — 1 communication cycles, but a hi- 
erarchical scheme (see table 3) sends the same amount of 
data with only ln2(n p ) communication cycles. The differ- 
ence between the two approaches remains small, as long as 
the communication is bandwidth limited, i. e. the blocks are 
large. Small block sizes shift the bottleneck to the latency, 
which is significantly reduced by the second scheme - espe- 
cially if the code runs on many processors. 



Table 4. Timings on a Beowulf cluster (Hydra, sec table 5). See 
text for an explanation of the variables. Timings are obtained 
for a maximal neighbour number LMAX=64. In practice, B is 
twice as large due to storage rearrangements in Nbody6++. See 
Appendix 5 for details on the computers. 
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64 
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64 


0.46 


401 
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Table 6. Timings on the IBM. More than 32 processors require 
more than one node. 



A hierarchical scheme reduces the latency, but never- 
theless it is possible that the parallel integration is actually 
slower than a single CPU integration. We estimated both 
the runtime on one CPU and on a parallel machine to ex- 
plore the transition between these two regimes. The latency 
time T( per communication is included in the wallclock time 
expressions for one regular/irregular step: 



T; = a A 

^single = aiVbl-Nnb 

A^bliVnb 



+ A\n 2 (n p )+ BNu 

Latency Communication 



(102) 
(103) 



If tsingic (runtime on a single CPU) is equal to t paT (parallel 
computation), one can deduce the critical block size N m - ln 
which gives the minimal block size for efficient parallelisa- 
tion: 



^single 
AUin 



(104) 



(-par 

Aln2(n p )n p 
N nh (n p — 1) — Bn p 

The hierarchical communication gives a minimal block size 
that increases logarithmically with the processor number. 
Eq. 103 gives immediately the speedup S and the optimal 
processor number for a certain block size JVbi: 



1 4- An 1 " 2(n P ) -1- R "P 



ln(2)JV M x AU 
n P! o P t(Abi) = -. — Hierarchical 



(105) 



A comparison to the optimal processor number for a ring 
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Name 


Beowulf (Hydra) 


Titan 


JUMP 


Institute 


ARI/ZAH 


ARI/ZAH 


Forschungszentrum Jiilich 


Location 


Heidelberg 


Heidelberg 


Jiilich 


Processors 


20 


64 


1248 


Speed 


2.2GHz 


3.2GHz 


1.7GHz 


Processors /Node 


2 


2 


32 


Network 


Myrinet 


Infiniband 


Gigabit— Ethernet 


Bandwidth 


2Gbit/sec 


20 Gbit/sec 


10 Gbit/sec 



Table 5. Specs of the different supercomputers used in running the algorithm 
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(106) 



stresses the efficiency of the hierarchical communication, 
since it allows a much larger processor number for a given 
problem size. Equation 102 and 103 are also useful to derive 
the total wallclock time, since the total runtime scales with 
the number of regular and irregular blocks: 



ATi/3 




(107) 



(108) 



These equations are only approximate expressions, but they 
give the right order of magnitude without detailed calcula- 
tions that need a precise knowledge of the A-body model. 
Table 4 and table 6 summarise the timing parameters drawn 
from our experience with the Hydra and JUMP parallel su- 
percomputers. 



8.6 Block Size Distribution 

The preceding section showed that the block size is closely 
related to the efficiency of the parallelisation. Small blocks 
are dominated by the latency and the parallelisation could 
be even slower than a single CPU calculation. Therefore we 
derive the block size distribution for the block time step 
scheme to asses its influence on the efficiency. 
Suppose that the time steps 14 h of all N particles in the 
model are distributed according to some known function /: 



dp = f(N, h)dh 



(109) 



/ is in most cases a complicated function. It involves spa- 
tial averaging and integration over the velocity distribution, 
which could be quite complicated even for simple time step 
formulae. Nevertheless there is a constraint on the time step 
distribution, simply because every particle has a neighbour 
within a finite distance: There is some upper limit /i max that 



dpi 

dh 



max 




Figure 5. Timestep distribution / = dp/dh. The short-dashed 
line on the left indicates approximation Eq. 110, whereas the 
dashed line on the right defines a reasonable upper limit /i ma x- 



restricts the major fraction of the time steps to a finite in- 
terval. Thus it is possible to capture the main features of 
the time step distribution with an expansion around h = 
(Fig. 5 sketches this approximation): 



/ « C(N)h a h^h n 



(110) 



a is the lowest non-vanishing order of the expansion. Now 
we consider a block level with the largest possible time step 
hk- The number of particles A^i in this block is: 



Abi = N / fdh i 



C(N) 
a + 1 



{h k 



(111) 



According to the block time step scheme the number of 
blocks per time with the largest possible time step hk is pro- 
portional to (hk) -1 . Therefore the probability that a block 
size is in the range [Nhi , Nu + dNbi] is 



dp<xJ2s{N hl -^^h a+1 



a + 1 



-dNu 



(112) 



14 We use h instead of At in this section to avoid unclear nota- 
tion. 



where <5 is Dirac's delta function. The sum over the loga- 
rithmically equidistant time steps hk is approximated by an 



Hybrid methods: A new composite algorithm 17 



3 



o.i 



0.01 



0.001 



le-04 



Plummet Model N=2500 
Theory a=2 




le-04 



0.001 



0.01 

N/N,. 



0.1 



Figure 6. Cumulative irregular block size distribution for a N ■ 
2500 particle Plummer model. 



integral: 



a + l 



(113) 



Thus the average block size and the median of the block size 
distribution are: 



(TVbi) « iiV o/(o+1) 



median (A^b 



(114) 



Special expressions for the average block size were already 
derived by Makino (1988), but the general relation of the 
time step distribution to the block size distribution is a new 
result. The median is surprisingly independent of the parti- 
cle number, i. e. 50 % of all blocks are always smaller than 
a fixed value. It seems that this is a threat to the efficiency 
of the method, but the median of the wallclock time 



median(iV bl ) 



N 



(115) 



demonstrates that these small blocks account only for a 
small fraction of the total CPU time. We confirmed the de- 
rived block size distribution (Eq. 113) by numerical calcula- 
tions (see Fig. 6). The order parameter a is roughly two in 
(at least locally) homogenous systems, while an additional 
Keplerian potential reduces the order to a = 1. A planetesi- 
mal disc - or more precisely, a narrow ring of planetesimals 
- has a very narrow distribution of time steps since all par- 
ticles share nearly the same orbital period. Thus the regular 
block size is always equal to the total particle number mak- 
ing the parallelisation very efficient. 



8.7 Optimal Neighbour Number 

We treated the mean neighbour number iV n b so far as some 
fixed value. But it is also a mean to optimise the speed of 
the integration. Large neighbour spheres reduce fluctuations 
in the regular forces allowing larger regular steps, which 
reduces the total number of force evaluations. But larger 
neighbour lists also imply a larger communication overhead, 
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Figure 7. Optimal neighbour number as function of particle 
number N. The plot includes the numerical solution of Eq. 124 
and the two asymptotic solutions. Timing constants are taken 
from a Beowulf cluster. 



as all the neighbour lists have to be broadcast to synchro- 
nise the different nodes. The best choice balances these two 
extremes, thus maximising the speed. 

Before we derive the optimal neighbour number on a 
parallel machine, we briefly summarise the known solution 
for a single CPU run (Makino 1988) for an extensive deriva- 
tion. The computational effort of the irregular steps is pro- 
portional to the neighbour number, while the number of 
force evaluations for the regular steps is proportional to the 
total number of particles, reduced by the time step ratio 7*: 



7t := 



At,, 



Afc, 



Tcpu = /(AO ( iVnb + 
7*(JV»b) « N^ 3 



N 



7t(N n 



(116) 
(117) 



f(N) collects all factors depending only on the total num- 
ber of particles. Optimisation with respect to the neighbour 
number iV n b yields the well known result: 



= 



d 



rib, opt 



dA/nb 

ociV 3/4 



Tcpu 



(118) 



The calculation of the elapsed time for NbodyGH — h on a 
PC cluster includes more terms. For clearness, we restrict 
ourselves to a rather simple model that involves only the 
dominant terms to show how parallelisation influences the 
optimal neighbour number. We make the following approx- 
imations: 

(i) We only take the force calculation and communication 
into account. 

(ii) We use the same time constants for regular and irreg- 
ular expressions. 

(iii) We neglect all numerical factors that are comparable 
to unity. 

The total CPU time is an extension of Eq. 106, which is 
applied to the regular and the irregular step. A new constant 
B n includes the neighbour list communication separately, 
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while all factors depending on N are represented by f(N): 
N bl « N 2/3 

T irr = /(iV) ( +Ap + BN hl 

T rog = -/(JV) + + (B + B„A/ nb )A/ b i 

7* V P 

Ttot = T irr + T rog (119) 

Optimisation with respect to the processor number p leads 
to: 

d m 
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TV 5 / 6 



(120) 



Further optimisation with respect to the neighbour number 
gives the expression: 

d 
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= A# » 



~-N- 



AN~ 2/3 p 2 



;Bp + 
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(121) 



(122) 



For a fixed p or B n = (very fast neighbour list communi- 
cation), we recover for large N: 



iVnb.opt OC N 3/4 



(123) 



In general, one can not neglect the neighbour list commu- 
nication. Therefore we seek for the optimal choice of p and 
7V n b, thus combining Eq. 120 and Eq. 122: 



ANT + ( jjiVn, 



VjNA 



(124) 



Since this equation has no closed solution, we identify 
the dominant terms in Eq. 124 to calculate the asymptotic 
solution for large N: 



iVnb, 



opt 



A 



3/5 



N 



1/5 



N>> {m 



3/2 



For small N we get the approximated solution: 

. ,3/5 / -1A \ 3/2 

■?Vnb,opt 



,r f 3AY 

N< {m) 



(125) 



(126) 



Fig. 7 compares the approximate expressions with the nu- 
merical solution of equation 124. In spite of the complicated 
structure of Eq. 124, both approximate expressions are re- 
liable solutions. The example uses timing constants derived 
form our local Beowulf cluster: 



200 B w 5 B„ 



0.5 



(127) 



If we compare the new optimal neighbour number to the 
single CPU expression (Eq. 118), we find that the influence 



of the neighbour list communication favours much smaller 
neighbour numbers. A^ n b increases so slowly with the particle 
number that a neighbour number around 100 is a safe choice. 



9 COLLISIONAL AND FRAGMENTATION 
MODEL 

The growth of planetesimals proceeds through collisions 
among planetesimals which form (at least in a sufficient frac- 
tion of incidents) larger bodies with a net gain of accreted 
matter. But some collisions are mere destructive events that 
shatter and disperse the colliding planetesimals. Small bod- 
ies are more susceptible to destruction, but they are also 
driven to high relative velocities due to the global energy 
cquipartition making them even more vulnerable. A model 
that attempts to cover the full size range from objects rang- 
ing between a kilometer and the size of Mars needs a realistic 
collision algorithm that covers both fragmentation and ac- 
cretion. Some examples in the literature are Cazenave et al. 
(1982); Beauge & Aarseth (1990). In our algorithm we use 
the approach of Glaschke (2003), which was applied to as- 
teroid families. The following section describes the approach 
to fragmentation that is implemented in the code. 

9.1 Handy quantities for quantifying the models 

Two colliding bodies are equal in the sense that their intrin- 
sic properties are not different. Only the comparison of two 
bodies defines the larger body - usually denoted as target 
- and the smaller one denoted as projectile. The two terms 
stem from laboratory experiments where they indicate much 
more than different sizes. A small projectile is shot on a 
target at rest to study the various parameters related to 
fragmentation. In the following, projectile and target only 
indicate the relative size of the two bodies. 

The collision of two bodies initiates a sequence of com- 
plex phenomena. Shock waves run through the material, 
flaws start to grow rapidly breaking the bodies apart in 
many pieces. Some kinetic energy is transferred to the frag- 
ments, which leads to the ejection of fragments at differ- 
ent velocities in various directions. If the fragment cloud is 
massive enough, some of the larger fragments may capture 
debris. This post-collisional accretion is denoted as reaccu- 
mulation. 

Although the depicted scenario is quite complex, there 
are a few measures that describe the most important as- 
pects: 

(i) Mass of the largest fragment Ml, or dimensionless 
fi = Ml/M where M is the combined mass of the two 
colliding bodies. 

(ii) fi < | refers to fragmentation, whereas /; > \ is 
denoted as cratering. 

(iii) Energy per volume S that yields /; 
as impact strength. 

(iv) /ke := 2_E^ s /_E kin : Fraction of the impact energy 
that is converted into kinetic energy of the fragments. 



| is denoted 



Different fragment sizes and velocities are summarised by 
appropriate distribution functions, m;, Di and Vi are mass, 
diameter and modulus of the velocity of a given fragment, 
respectively. 
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Figure 8. Section for /; = 0.04 and n = 3. The largest fragment 
is coloured in dark-grey. In this calculation 60 X 60 X 60 grid cells 
are used. Note the decomposition in grid cells and the Voronoy 
polyhedra which form the fragments. 



We selected two different impact strength models as 
reference for our work. The first was obtained by Housen 
(1990) through the combination of asteroid family data and 
laboratory experiments via scaling laws: 



(i) Fragment size distribution: 

(a) N m (m) : Number of all fragments with a mass rrt; ^ 

m, 

(b) M(m) : Mass of all fragments with a mass mi ^ m, 

(c) Nd{D) '■ Number of all fragments with a diameter 
D, > D. 



The distribution functions are related to each other: 
N m (m) 

M(m) - I r ' ;'V' J ,/.,. 



N D (D(m)) 

dN m (x) 



N m {m) = 



dx 
dM(x) 



dx 



dx 



D(m) is the size-mass relation, 
(ii) Velocity distribution: 

(a) v(m): mean velocity as a function of mass. 



9.2 Prediction of collisional outcomes: Derivation 
from a Voronoy tessellation 

Any theoretical or empirical prescription of a collision has to 
relate the afore mentioned parameters, namely the impact 
energy, to the sizes and velocities of the produced fragments. 
The central quantity is the impact strength, which is a mea- 
sure for the overall stability of a body. Objects smaller than 
1 m are accessible to laboratory experiments, while collisions 
of larger bodies up to asteroid size have to be analysed by 
complex computer simulations. Asteroid families, which are 
remnants of giant collisions in the asteroid belt, provide in- 
dependent insight, although the data is difficult to interpret. 



So 



/kb = 0.1 



R 
lm 



1 + 1.6612 x 10" 7 



R 
lm 



So = 1.726 x 10 6 Jm" 3 = 1.726 x 10 7 erg cm" 



(128) 
(129) 



Later, Benz (1999) obtained another result through 
SPH simulations (for basalt, v —3 km/s): 



S 

/ke ; 



So 
0.01 



R 
lm 



1 + 6.989 x 10" 



R 
lm 



So = 6.082 x 10 s Jm" 3 



2.7- 



(130) 
(131) 
(132) 



/kb is a measure of the kinetic energy that is transferred to 
the fragments: 



frag 



"-Ekin 



(133) 



We introduce a dimensionless measure 7 of the relative im- 
portance of gravity for the result of a collision. It is defined 
as the ratio of the energy per volume Sg that is necessary 
to disperse the fragments to the impact strength So: 



S G = 2tt- , 

5/ke 

7 := Sg/S 



2 ^GR 2 p 2 



(134) 



The first step towards the prediction of a collisional out- 
come is to relate the impact energy and the impact strength 
to ascertain the size of the largest fragment Laboratory 
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experiments and simulations indicate the functional form 



<fi) = 



2(1-/0 for /,>§ 
(2fi)~K otherwise. 

2SM 



(135) 



which is both valid in the fragmentation regime and the cra- 
tering limit. The size of the largest fragment is used to derive 
the full size distribution. To accomplish the decomposition 
"seed fragments" are distributed inside the target accord- 
ing to the largest desired fragment. The full set of fragment 
is derived from a Voronoy tessellation 15 using these seed 
points. Fig. 8 depicts the result of such a decomposition. 
The fragment velocities are calculated from the total kinetic 
energy after the collision to initiate a post-collisional N— 
body calculation to treat reaccumulation. 

We conducted a large set of such calculations to cover a 
sufficient range in / ( ' (i. e. impact energy) and 7 (i. e. body 
size). Table 7 summarises the derived values of the largest 
and second largest fragment including reaccumulation. 



10 COLLISIONAL CASCADES 

A first well-defined application of the fragmentation model 
is a collisional cascade. The term cascade denotes that frag- 
ments of one collision in a many-body system may hit other 
bodies, whose fragments further shatter even more bodies. 
Thus the particle number increases exponentially with every 
subsequent collision. 

Although the formation of planets requires a net growth 
due to collisions, this destructive process plays a role in the 
formation of larger bodies as the overall size distribution 
controls the accretion rate of the protoplanets. Therefore it 
is worth to have a closer look into this mechanism. 



10.1 Self— similar collisions 

A system of colliding bodies is usually embedded in a 
broader context, like stars moving in a galaxy or asteroids 
orbiting in our own solar system. First, we simplify this dy- 
namical background as well as some aspects of the collisions 
to make the problem tractable. 

The first step is to decompose an inhomogeneous sys- 
tem into smaller subvolumes which are locally homogenous. 
Furthermore, it is assumed that these subvolumes hardly 
interact with each other. Hence it is possible to apply the 
particle-in-a-box-method (Safronov 1969) to analyse colli- 
sions within the small subvolumes: 

(i) All particles are contained in a constant volume. 

(ii) The particle sizes are described by a distribution func- 
tion n(m), i.e. the particle number per volume and mass 
interval. 

(iii) For convenience, we assume a constant (or typical) 
relative velocity for a given pair of colliding bodies. 



15 The Voronoy tessellation assigns every volume element to the 
closest seed point. First applications date back to the 17th cen- 
tury, but the Russian mathematician of Ukrainian origin Gcorgy 
Fcodoscvich Voronoy put it on a general base in 1908. 



The distribution function is evolved by the coagulation equa- 
tion. We modified the equation given by Tanaka et al. (1996) 
by introducing a new function M re d to arrive at a more con- 
cise expression: 



d d 
= — mn(t, m) + ——F m (t, m) 
dt dm 



(136) 



The mass flux F m is given by: 



M, 



m = — J J n(t,m\)n(t,rn2)(,dm,idm,2 
£ = <7(mi,m 2 )«reiM re d(m, mi,m 2 ) 
ot = J n(t, m)mdm 



d_ 

dt 



Mtot — F„ 



(137) 
(138) 

(139) 



where £ is the coagulation kernel, n is the already intro- 
duced size distribution, v Te \ is the mean relative velocity, a is 
the cross section for colliding bodies (mi, 7712) and M ro( j is 
the newly introduced fragment redistribution function. M rc d 
contains all information on the fragments arising from the 
breakup of body mi due to the impact of body m 2 . Its defi- 
nition avoids double counting of collisions in the above inte- 
gral. The redistribution function is related to the differential 
number distribution function n co ii(mi, m 2 , m), i. e. the num- 
ber of fragments produced by a collision per mass interval. 
Since the target mi formally disappears, it is included as a 
negative contribution: 



M red (m,mi,m 2 ) := 

f-ra 

I (ricoii(mi,m2,m) — 5(m — m{))mdrh 
Jo 



(140) 



Mass conservation in each collision is reflected by 
M r ed(0, mi, TH2) = M re d (00, mi, m-i) = 0. The cross section 
a depends on the velocities and radii R t of the particles. A 
simple approach is the geometric cross section: 



a(mi, m 2 ) — tv(Ri + R2) 



(141) 



If gravity plays an important role during encounters, two 
colliding bodies move on hyperbolic orbits with a pericentre 
distance that is smaller than the impact parameter. This 
leads to an additional enlargement of the cross section, de- 
noted as gravitational focusing: 



a(mi, m 2 ) = n(Ri + R2) 2 ( 1 + 



2G(mi 



m 2 , 



l (Ri+R 2 ) 



(142) 



A special class of collisional models are self-similar collisions. 
Self-similarity implies an invariance of the collisional out- 
come with respect to the scale of the colliding bodies. If the 
target mass as well as the projectile mass are enlarged by a 
factor of two, then only the masses of all fragments doubles 
without further changes in the collisional outcome. They al- 
low us to introduce the convenient dimensionless fragment 
redistribution function f m : 

M rcd (m, mi,m 2 ) = mf m (m 1 /m, m 2 /m) (143) 

We follow Tanaka et al. (1996) and employ the substitu- 
tion 16 mi = mil, m 2 = mx2 to simplify Eq. 137: 



16 A similar approach to the solution of the coagulation equation 
is the Zakharov transformation, see Connaughton et al. (2004). 
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Largest Fragment 



7 ft 


0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


0.8 


0.02552 


0.10000 


0.20000 


0.30000 


0.40000 


0.50000 


0.60000 


0.70000 


0.80000 


0.19897 


0.10000 


0.20000 


0.30000 


0.40000 


0.50000 


0.60000 


0.70000 


0.80000 


0.67985 


0.10000 


0.20000 


0.30000 


0.40000 


0.57612 


0.67315 


0.82006 


0.96073 


1.14050 


0.10000 


0.31526 


0.35708 


0.61362 


0.83511 


0.92832 


0.94380 


0.97572 


1.77057 


0.10884 


0.58883 


0.75974 


0.87922 


0.92755 


0.93662 


0.97107 


0.97924 


2.26021 


0.15895 


0.68891 


0.87217 


0.89592 


0.92965 


0.96089 


0.96701 


0.98727 


3.11626 


0.30954 


0.83774 


0.90272 


0.92682 


0.95943 


0.95565 


0.97677 


0.98791 



Second largest Fragment 



7 ft 


0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


0.8 


0.02552 


0.08171 


0.10770 


0.06471 


0.08107 


0.04976 


0.03982 


0.03171 


0.02155 


0.19897 


0.09174 


0.09533 


0.08410 


0.06967 


0.04930 


0.04825 


0.03510 


0.02077 


0.67985 


0.07713 


0.08365 


0.07791 


0.08387 


0.07026 


0.06147 


0.06082 


0.00847 


1.14050 


0.08621 


0.07256 


0.10331 


0.09640 


0.02467 


0.00675 


0.00783 


0.00265 


1.77057 


0.07909 


0.05549 


0.06961 


0.02035 


0.00329 


0.00719 


0.00273 


0.00161 


2.26021 


0.06693 


0.02288 


0.00528 


0.00882 


0.00584 


0.00268 


0.00664 


0.00126 


3.11626 


0.06940 


0.00884 


0.00384 


0.00488 


0.00064 


0.01007 


0.00225 


0.00162 



Table 7. Data compilation of the fragmentation calculations. 



to Eq. 136 using the ansatz n(t,m) = a(t)no(m): 



F m = — J j n(t,rnxi)n(t,m,X2)m 11 ^ i 

a(xi,x 2 )v lcl f m (xi,x 2 )dxidx2 (144) 

A simple solution is a steady-state cascade with F m — const. 
The loss of bodies of a given size is balanced by the frag- 
ment supply from larger bodies, hence the system maintains 
a steady-state ^n(t,m) = 0. Eq. 144 inspires the ansatz 
n(m) oc mT k , which yields k = 11/6. This is the well known 
equilibrium slope in self-similar collisional cascades, which 
was already found by Dohnanyi (1969). Strong gravitational 
focusing changes the exponent 17 to k — 13/6. Both steady- 
state solutions seem to be rather artificial, as they contain 
an infinite amount of mass and require a steady mass influx 
from infinity. However, they provide an appropriate descrip- 
tion for the relaxed fragment tail of a size distribution, as 
long as the largest bodies provide a sufficient flux of new 
fragments. Once the largest bodies start to decay, the finite 
amount of mass in the system leads to an overall decay of the 
collisional cascade. Thus we seek for a more general solution 



17 Tanaka et al. (1996) state that k < 2 is a necessary condition 
for a finite mass flux. However, their analysis is not valid for all 
possible collisional models. 



-Ca(ty 



mn (m) = -^-^—F„, 
C dm 



(145) 
(146) 



C is determined by fixing no at an arbitrary value m* . a(t) 
is independent of the collision model: 



a(t) 



l + Ct 

C oc n(m*) 



(147) 
(148) 



A power law solution is no(m) oc -Cm _fc+1 which is only 
valid for C < (agglomeration dominates). To examine C > 
0, we perturb the already known equilibrium solution: 



no(m) 
1 

ivT 



N m~ k 
(2-fc) 



-2k+2 
-2k+2 



a(xi,X2) 



Vrel {fm(xi,X 2 ) + f m (x2,Xl)) dxidx 2 



(149) 



(150) 



Ni is small if the integral on the right hand side is large. This 
is the case for a sufficiently large impact strength. Eq. 149 
has the interesting property that n(m') = for some mass 
m', given that k < 2. This mass m' represents the largest 
body in the system, e. g. the largest asteroid in a fictitious 
asteroid belt. 
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11 SIZE-DEPENDENT STRENGTH 

Self-similarity is an enormous help in analysing the coagu- 
lation equation. It releases completely the need to know any 
specific details of the collisional process and provides valu- 
able insight at the same time. But self-similarity is also a 
strong limitation on the underlying collisional physics. 

A major component of a fragmentation model is the 
knowledge of the impact strength as a function of size. Sim- 
ulations as well as asteroid families establish that it is not 
some fixed value, but changes with size which immediately 
breaks the self-similarity. Larger bodies are weaker due to 
an increasing number of flaws (there are no big monocrys- 
tals), but then gravity leads to a turnover and increases the 
strength. 

We model the size dependent strength S with a power 
law to examine the influence on the equilibrium solution. 
The velocity dispersion v and the collisional cross section a 
are also modelled with power laws to account for relaxation 
processes: 



(-)' 

m \ ' 



V = vo 



a = cro I — 



S = S 



(151) 
(152) 
(153) 



The subscript "0" denotes values for an arbitrarily chosen 
scaling mass. Since smaller bodies are more abundant than 
larger ones, we safely assume that most collisions involve a 
large mass ratio. In addition, we assume w < 0, since we 
expect energy equipartition to some degree in most cases. 
These restrictions lead to the following simplifications (mi > 
m 2 ): 



cr(mi, 7712) ~ a(roi) 
?j ro i « v{m 2 ) 

~ 2 mi Si 



(154) 
(155) 

(156) 



Therefore the smaller body WI2 enters only through the spe- 
cific energy e: 

F m « - J J n(t, mi)n(t, m2)a(m 1 )v Icl (m 2 )m 1 

f m (mi/m,e)dmidm2 (157) 

We introduce new dimensionless quantities with the help of 
Eq. 154-156 to simplify the integral: 



mi = mxi 



m,2 = mo 



/mi \ 1 
\m J 



1 1 ■< 
. 2. 



/ 2So \ _l_ 
— £ e'+2» 

\PvlJ 



(158) 



Again we assume a power law for the density n oc m~ and 
change the integration parameters to (xi,e). Applying the 
constant -flux condition yields the equilibrium exponent 



k 



S + 3 + Q + w(2s + a + 5) 
2 + a + 2w 



(159) 



and the scaling exponent k' of the total mass loss: 
s — w + 1 



k' 



dt 



tot OC 



2 + a + 2w 

-k' ? 



2So 

PVq 



(160) 



The exponent k' in Eq. 160 is close to unity for realistic 
values of the free parameters. Thus the mass loss is roughly 
inversely proportional to the strength of the bodies. The gen- 
eral formula Eq. 159 contains the special solution of O'Brien 
(2003), who concentrated on the parameters s = 2/3, w — 
and a special collisional model. In fact, the derivation applies 
to a much wider class of collisional models that we denote as 
scalable collisional models. Scalable indicates that the model 
is self-similar except a scaling of the impactor mass. 



12 PERTURBATION OF EQUILIBRIUM 

The derived scaling relations provide insight into the over- 
all properties of a collisional cascade, which is in (or close 
to) equilibrium. However, they do not provide information 
on how the equilibrium is attained or how the system re- 
sponds to various external perturbations. A rigorous ap- 
proach would be the approximate solution of the coagulation 
equation 18 , which is by no means simple since it requires a 
careful analysis of the collision model. 

Hence we turn to perturbations of the equilibrium size 
distribution, as it is easier to asses the quality of the derived 
expression for a variety of collision models. In addition, all 
equations are linear in the perturbation, allowing the de- 
tailed analysis of the solution. 

If the equilibrium solution n(m) = no(m/mo)~ k is per- 
turbed with a small deviation An(m), we get to first order: 



0: 



d d 
— mAn(m) + —F p (t,m) 



(161) 



J I An(mi)n(t, m 2 )a(m 1 , m 2 )v Icl 



x (M rc d(m, mi, WI2) + Mred (m, m2, mi)) dm\dm2 

Despite of the expansion in An, Eq. 161 is still a compli- 
cated integro-differential equation. Thus it is not possible 
to obtain a solution without further information about the 
problem. While there is no general solution, we restrict our 
attention to self-similar collisional processes. In virtue of this 
assumption it is possible to simplify Eq. 161, as we can see in 
Eq. 162 and 162. In those expressions cro and vo are velocity 
and cross section of an arbitrarily chosen scaling mass mo- 
-F(:ei) contains all information about the collisional process. 
If collisions do not result in extreme outcomes, like cratering 
or a complete destruction of the target, most of the fragment 
mass is contained in bodies with similar size as the parent 
body. Hence we expect that F(x\) peaks around xi ~ 1 and 
drops to zero as X\ gets larger (or smaller). We introduce 
the dimensionless relative perturbation g(m): 



An(m) _ An(m)m fc 
n(m) 



(164) 



Appendix C highlights a possible approach. 
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o = — mA 
dt 



n(m) — norrtQcrofo f An(t, mxi)F(xi)(mxi /mo)' dx\ 

dm J 



EV \ f 2k-3 -k -k a ( x l> x 2) V ie l , . 
F(xi) = / m x 1 x 2 —{fm{xi,x 2 ) + fm{x2,x 1 ))dx 2 



'■() 



(162) 
(163) 
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Figure 10. Scaled fragmentation kernel G(u) for a simple 
fragmentation model (see Eq. 200) and different scaled impact 
strength S. 



Thus the new differential equation reads: 

d 

ot 

2 d 
n m Q a v -^ 



Ag(t, mxi)F(xi)dxi 



(165) 



We change to logarithmic coordinates to arrive at a convo- 
lution integral: 

u — ln(m/mo) u\ = \n(x\) (166) 

Furthermore we define a collisional timescale to 

t = (nomoGQVoY 1 (167) 

to obtain a more concise expression. The transformed equa- 
tion is: 



g(t,u) 



u(2-k) 



" dt 



G(u) =F(e u )e J 



1JL 
to du 



g(t, u + ui)G(ui)dUi 



(168) 
(169) 



If g(u) is varying on a scale larger than the width of the 
kernel G(u) (compare Fig. 10), it is justified to expand g(u) 
under the integral. We retain the first two moments of G(u): 



Gi & 



n ^ \ u(2-k) Go 8 . . 

Cft To OU Tq OU 



a 2 f( t ' w ) 



G k 



u G(u)du 



(170) 
(171) 



The first order moment Gi , which introduces a diffusive 
term, is omitted in the following for clarity 19 . We introduce 



a fragmentation time Tf rag (M) and transform Eq. 170 back 
to m: 



= §i9(t,m) 



d 



Tf r ag(m) dm 



g(t,m) 



'frag 



TO u(2-k) 

(jo 

To 



Go 



(m/mo) 



(172) 



(173) 



Eq. 172 is a modified advection equation, which conserves 
the total mass. It is possible to derive equations similar to 
Eq. 172 for any collisional model. However, the general ap- 
proach is less fruitful, as it lacks a robust frame of a known 
equilibrium solution and reliable scaling relations. Therefore 
we provide only the extension to scalable collisional models 
in Appendix B. We readily obtain the general solution: 

(m/m ) (2 - fe) ' 



g(t, m) = f [t + T 



An(t, m) — n(m)f I t + To 



G (2-fc) 

(m/m ) (2 - fc) 



(174) 



G (2-fe) 

The function / is determined by the initial value g(0, m) 
of the perturbation. As the collisional cascade evolves, the 
initial perturbation function is shifted as a whole to smaller 
masses. This evolution becomes clearer if we attach labels 
M (0) to the initial perturbation function and follow the time 
evolution of these tags. The functions M(t) are the charac- 
teristics 20 of the differential equation 172: 



(t) = mo [(A/(0)/m ) (2 - fc) - */t G (2 - 



M 



The meaning of the fragmentation time Tfe 
by the relation 

M 
M 



Tf'rag 



-k)] 1/(2 - k) 

(175) 
becomes clear 



(176) 



which is the time until a body has lost a significant fraction 
of its mass due to destructive collisions. A comparison of the 
perturbation equation 172 with the scaling relations from 
the previous section gives the scaling of the zeroth order 
moment Go with respect to the impact strength: 



Go — G'qS 



(177) 



G'o should only depend on the fragmentation model (i. e. 
fragment size distribution as a function of the largest frag- 
ment /() within the limits of this approximation. Fig. 10 
shows that the scaling with the impact strength works quite 
well, except slight variations which are small compared to 
the covered range of impact strengths. Likewise, it is pos- 
sible to restate the total equilibrium flux F ecl in terms of 
G'o: 

~ Go 



F a 



-n(m) a(m)m v Te iS 



(178) 



19 The study of wave-like structures in the size distribution 
(Bagatin et al. 1994, see e.g.) requires even the second order mo- 
ment G2- 



20 In general, characteristics of a partial differential equation are 
paths along which the solution is constant. 
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The fragmentation timescale Tf rag (m) allows a more intuitive 
expression: 

^(m)«-i=^ (179) 

2 Tfrag (m J 

Our simple collisional model (see Fig. 10 and Eq. 200) refers 
to: 

Feq(m) = -(1 . . . 30) x n{m) 2 a{m)m 3 v Tcl S' k ' (180) 



13 MIGRATION AND COLLISIONS 

The local perturbation analysis is only applicable to a plan- 
etesimal disc, if the migration velocity of the planetesimals 
is negligible small. This assures that collisional cascades at 
different radial distances do not couple to each other, so that 
the whole disc is composed of many local cascades. While 
this assumption is justified for larger bodies, migration is 
strongly influencing bodies below 1 km in size. Hence we 
extend our analysis to examine the influence of migration 
on the (no longer) local collisional processes. 

We assume that the collisional evolution of the system 
leads to an equilibrium planetesimal distribution everywhere 
in the disc: 



E (r,m) = E ri0 (r)C (m) 



(181) 



E r (r) is the total surface density at a given distance r, while 
Co(m) is the universal equilibrium distribution. Though the 
planetesimal distribution at larger sizes is likely different at 
different locations in the disc, we only demand a universal 
function at smaller sizes, where migration is important. The 
power law exponent k depends on the details of the invoked 
physics, but numerical simulations show that k « 2 is a 
fiducial value. Eq. 181 does not yet include migration effects. 
If we include migration, the surface density is modified to 



E(r,m) = g(r,m)Y, (r,m) 



(182) 



where the dimensionless function g contains the changes due 
to migration. The collisional evolution is governed by the 
continuity equation with an additional collisional term 

<9E(r,m) 1 d 



(v(r, m)rT,(r, m)) 



(183) 



dt r dr 

where v(r, m) is the migration velocity (see Eq. 77), defined 
such that positive v imply an inward migration. We express 
the collisional term with the help of Eq. 172 and seek for a 
steady-state solution E = 0: 

—J r !^mE r , (r) + i ~ (pwrB r , (r)) = (184) 

Tfrag (m, r) am r or 

T£r a g(m, r) is the fragmentation timescale of a mass m at 
a distance r. Since the surface density E and the various 
contributions to the drag force are well described by a power 
law (with respect to radius), Eq. 184 further simplifies to: 

b 



1 dg ^ dg ^ 
Tfrag (m, r) dm dr 



gv 







(185) 



b is a combination of the various invoked power law 
exponents. As the surface density E and the gas density 
drop with increasing radius in any realistic disc model, it is 
safe to assume b > 0. We choose a self-similar ansatz for g: 




Figure 11. Cut-off function g according to Eq. 192. The mass 
exponent is k m = 1/3, while the mass influx exponent is b = 1.75 
according to the minimum solar nebula. 



The new differential equation is 

1 dq , . dq dq n 
mg m (r) +m-f--* 



v gv 







(187) 



(188) 



Tfrag(m,r) dC ' ' "~d( dr 

which is equivalent to the more concise expression 
dln(g) ( r dln(g m ) \ _ fe 
rfln(C) V wr frag dln(r) / 
We assume a power-law dependence for the timescale ratio 

Tmig / Tfrag - 

^_ = ZiE*. = (m/m ) fcm (r/r- ) fc " (189) 

^Tfrag Tf ra g 

The cut-off mass mo at a distance ro has a timescale ratio 
Tmig/Tfrag = 1, which defines a proper lower cut-off within 
this context. Hence the solution is: 
h _dln(g) 



9m(r) = 



9(0 



dln(C) 

(r/r ) fc " /fc 



C m + 



m 



1 + 



k t k 



-b/k T 



(190) 
(191) 
(192) 



Though the analytical solution Eq. 192 provides a complete 
description of the lower cut-off of the size distribution, it 
is more appropriate within the frame of this discussion to 
translate the equilibrium solution to an equilibrium mass 
loss due to migration: 



kr / krt 



. , bv bv 
6E 



(193) 



g(r,m)=g(0 , C = rng m (r) 



(186) 



Tmig ~t~ k r I k m Tf ra g 

An inspection of the timescale ratio shows that the mass 
exponent k m should be positive, whereas simple estimations 
of k r on the basis of the minimum mass solar nebula are 
somewhat inconclusive. The value of k r is so close to zero 
that any change in the assumed equilibrium slope or the 
impact strength scaling gives easily both positive and neg- 
ative values. Moreover, Eq. 193 requires a globally relaxed 
planetesimal disc, but the huge spread in the various in- 
volved timescales at different radii inhibits any significant 
relaxation in the early stages. 
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However, it is possible to gain valuable information from 
the two limiting cases k r > and k r < 0. Both values of 
k r give the proper limit g — > 1 at large masses, where the 
migration timescale is much larger than the fragmentation 
timescale and we recover the steady-state collisional cascade. 

A positive exponent k r reduces the effective mass loss 
due to migration, as fragments from the outer part of the 
disc replenish the local mass loss. Hence the fragmentation 
timescale controls the net loss of smaller planetesimals. In 
contrast, a negative exponent k r leads to a pronounced cut- 
off in the size distribution, since only larger planetesimals 
are replenished through inward migration. Though the mass 
loss rate is singular at some mass m' , this sharp cut-off is 
an artifact due to the perturbation approximation. 

Our analysis is subjected to several restrictions. We ap- 
plied the perturbation equation to values of g that exceed 
the limit for a safe application (i. e. g 96 1) of the perturba- 
tion expansion. Furthermore, the steady-state solution re- 
quires a global relaxation of the collisional processes, which 
is practically never obtained during the disc evolution. De- 
spite of these restrictions, we gained insight on a more qual- 
itative level. Numerical calculations indicate that the per- 
turbation approximation is inappropriate close to the lower 
cut-off of the size distribution. However, a comparison of dif- 
ferent exponents k r (see Fig. 11) attributes only a minor role 
to the replenishment of fragments due to inward migration. 
Only unrealistic small slopes b of the migrational mass influx 
would strengthen the importance of this process. Though 
temporally non-equilibrium phenomena are not ruled out 
by the previous derivation, their study would require the 
global simulation of the system. 



14 COAGULATION 

While most coagulation kernels are only restricted to a lim- 
ited analytical analysis (e. g. scaling relations) , there exist 
some special kernels that allow the closed solution of the 
coagulation equation. All rely on the assumption of perfect 
mergers, which allows the reformulation of the general equa- 
tion 136 to 



dt 



dn[m,t) 1 f m , , , , , 

— - / /C(m — m ,m )n(m — m ,t)n(m ,t)am 
2 Jo 

■n(m,t) / K(m,m')n(m ,i)dm 
Jo 



where K, is the coagulation kernel. One of these particular 
kernels was introduced by Safronov (1969): 



K(mi,m 2 ) = Ai(m 1 + m 2 ) 



(194) 



This coagulation kernel implies perfect mergers, where the 
coalescence rate of two particles mi and m 2 is assumed to be 
proportional to the sum of their masses. It seems that this 
is an artificial choice, devised to allow an analytic solution. 
However, the Safronov cross section provides an intermedi- 
ate case between a geometric cross section (a oc m 2 / 3 ) and 
strong gravitational focusing (a oc m 4//3 ). A special solution 
to the initial condition 



n(m, 0) = cxp (— ml mo) 
m 



p 


2,700 kg/m 3 


k 


1/6 


Model 


Gaussian Scatter 


Ike 


0.1 


K 


1.24 



Table 8. Main parameters of the collisional model. 



is the function (Ohtsuki et al. 1990, see e.g 
n 



(m,r) = — " 0g = cxp(-(2 - g)m/m )J"i(2m/mo\A - g) 
my/1 - g 

(196) 



g = exp(-r) r = A L pt 

/>oo 

p = I mn(m)dm = nomo 
Jo 



(197) 



where r is the dimensionless time and h is a modified Bessel 
function of the first kind. 



15 MODELS FOR M RED 

Though we already obtained insight into the nature of col- 
lisional cascades without a detailed specification of the co- 
agulation kernel, any detailed study of a collisional system 
requires the specification of a realistic collisional model. 

First, we restate the well-known perfect accretion 
model. While it is a gross oversimplification for collisions 
among kilometre-sized planetesimals, its simplicity allows a 
reliable code testing and eases the comparison with other 
works: 

M le d(m, mi, m-i) = — miO(m — mi) — m20(m — m2) + 
(mi + m,2)6(m — mi — 7712) (198) 

Although our fragmentation model (see section 9) pro- 
vides a very detailed description of the outcome of a colli- 
sion, we abandon most of the details for the following rea- 
sons. The computational effort of the numerical solution of 
the coagulation equation scales with the third power of the 
number of mass bins. Hence we chose a mass grid whose res- 
olution is by far smaller than the information provided by 
the detailed collisional model. As a mismatch of the mass 
resolution could produce undesired artifacts, a lower resolu- 
tion of the collisional model is needed for consistence. Thus 
only the largest fragment fi(f{, 7) and the second fragment 
/j*- 2 ' (f{, 7) (which contains information on reaccumulation) 
enter the fragment size distribution: 

'1 ifx>/« 
1 



M rcd {xM) 



M 



fi 



(l-/i-/ ; (2) )(x// ( w ) /! otherwise 



if fi > x > // 2) 



(199) 



1,(2) — 

Both values /; and /; are interpolated from table 7, where 
the initial fragment size ft is calculated from the dimen- 
sionless impact energy e. We used a reduced fragmentation 
model for test purposes: 

MM*M)/M = { ;f^ //i)/!otherwise (200) 



(195) 
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Table 8 summarises the most important model parameters. 



16 STATISTICAL MODEL 

The direct approach to the integration of an iV-body system 
is, in principle, possible for any particle number. While this 
procedure becomes computationally too expensive for very 
large particle numbers, a by far more efficient approach is 
applicable in this regime. Instead of tracking all particle or- 
bits, a distribution function / (also phase-space density ), 
which gives the probability to find a particle at a position x 
with a velocity v, contains the state of the system: 



dp = /(x, v)d 3 xd 3 v 



(201) 



As long as only dynamical interactions are taken into ac- 
count, the number of all particles (e. g. stars, planetesimals) 
is conserved. The continuity equation reads: 



0-|£ + T.*/-V..fi 



(202) 



g + ,.v/-v..« 



(203) 



This is the collisionless Boltzmann equation. Collisions lead 
to an additional term 

'dp 

. 9t / coll 

which will be discussed later. / is a function of six variables, 
so an exact solution is usually very complicated or even im- 
possible. However, it is possible to gain valuable insight into 
the problem by taking the moments of the distribution func- 
tion: 



{XiVj ) = / / (x, vJrEj Vj d xd v 



n,m> (204) 



The spatial density (particles per volume) is related to dis- 
tribution function: 



i/(x) = j /(x,v)d 3 « 



(205) 



Integration of Eq. 203 over all velocities yields the corre- 
sponding continuity equation: 



dv dvvi 
dt dxi 



(206) 



The first order moment with respect to velocity gives the 
time evolution of the mean velocity v 



dvj 



dvj 



dt ~*~ VV% dxi 



dxj 



dxi 



+ v 



dt 



coll 



Vi = i J f{-x.,w)vid 3 v 



■ ViVj — ViVj 



(207) 



(208) 



where Oij is the anisotropic velocity dispersion and the con- 
tinuity equation was used to arrive at a more concise for- 
mulation. Equations 206 and 207 are the Jeans equations. 
While the structure of the moment equations is already fa- 
miliar from hydrodynamics, they do not provide a closed set 
of differential equations, since each differential equation of 
a given moment is related to (yet unknown) higher order 
moments. Hence any finite set of momenta needs a closure 
relation - additional constraints that relate the highest order 
moments to known quantities. The choice of this relation is 
a key element in the validity of the equations, but it is not 



unique and depends well on the problem at hand (Larson 
1970, compare e.g.). 

Owing to the geometry of a planetesimal disc, it is useful 
to express the Boltzmann equation in cylindrical coordinates 



dt r dr z dz \ r dr ) dv r 



VrV4, df _ gg df 
r dvs dz dv. 







(209) 



where all derivatives with respect to <f> have been dropped 
due to the assumed axisymmetry of the disc. 



16.1 Distribution Function 

Any statistical description of a planetesimal disc requires 
the knowledge of the distribution function. Since the full 
problem including collisions, encounters and gas drag has 
no analytic solution, a collisionless planetesimal disc (i. e. no 
perturbations) is a natural basis for further investigations. 
The distribution function that describes such a simplified 
system is a solution of the Boltzmann equation. A special 
solution to Eq. 209 is a thin homogenous planetesimal disc 



/(*,«) 



QE 

2-K 2 T r T z m 



exp 



2T r 



vl + QV 
2T 2 



(210) 

provided that the radial velocity dispersion T r and the ver- 
tical dispersion T z are small compared the mean orbital ve- 
locity vk- The azimuthal velocity dispersion T$ is locked to 
T r by the local epicyclic frequency k in a central potential, 
where the ratio 1 : 4 is a special solution of (Binney 1994, see 
e-g) 



(211) 



All velocities v r , v^, and v z refer to the local Keplerian ve- 
locity. The normalisation is the same as in Stewart (2000): 



vdzf(z, v) 



E 
m 



(212) 



A planetesimal disc is a slowly evolving system compared to 
the orbital time, hence it is reasonable to use Eq. 210 as a 
general solution of the perturbed problem. E, T z and T r are 
now functions of time and of the radial distance to the star. 
All information on the system is contained in these three 
momenta of the distribution function, where higher order 
moments can be deduced from Eq. 210. Thus the functional 
form of the distribution function represents an implicit clo- 
sure relation. 

The validity of this approximation can be further as- 
sessed by a closer examination of the Boltzmann equation. 
We summarise all perturbations in an evolution timescale 
T CV oi and reduce the radial structure to some typical length 
scale Ar to estimate the deviation from the functional form 
Eq. 210. A comparison with Eq. 209 shows that the differ- 
ence is small if the migration timescale and the evolution 
timescale are large compared to the orbital time T : 



To < Ar/(v r ) 

To <C T evo l 



(213) 
(214) 
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An order-of-magnitude estimate of the evolution time sup- 
ports condition 213 and 214. Furthermore, numerical calcu- 
lations confirm that the velocity distribution stays triaxial 
Gaussian (see Ida 1992) . 

The distribution function is equivalent to an isothermal 
vertical density structure with scale height h: 



h = 



p{z) = po exp 



z 

2h? 



(215) 
(216) 



Thus the central density p and the mean density (p) are 
related to the surface density in a simple way: 



Po 



(P) = 



V27T/1 

PQ 
V2 



(217) 



The triaxial Gaussian velocity distribution is equivalent to 
a Rayleigh distribution of the orbital elements e and i 21 : 



dn(e 2 ,i 2 ) = 



<e 2 > = 



{ei){i*) eXP { {*) 



de 2 di 2 



2T r 
WoT 2 



(i 2 ) 



2T Z 
(fin.) 2 



(218) 



Planetesimal encounters couple the time evolution of eccen- 
tricity and inclination, so that the ratio i 2 /e 2 tends to an 
equilibrium value after a few relaxation times. It is close 
to 1/4 in a Keplerian potential, but the precise value also 
depends on the potential itself (Ida et al. 1993). 



16.2 Dynamical Friction 

Planetesimal-planetesimal scatterings change the velocity 
distribution through two different processes. Firstly, it is 
unlikely that two planetesimals scatter each other on circu- 
lar orbits. Thus we expect a steady increase of the velocity 
dispersion due to this viscous stirring. Secondly, encounters 
between unequal masses lead successively to energy equipar- 
tition, slowing down the larger bodies through dynamical 
friction. The later mechanism is not related to the disc ge- 
ometry at all, but operates in any multi-mass system. A 
special case is the systematic deceleration of a massive body 
M in a homogeneous sea of lighter particles m with density 
no, which is given by the Chandrasekhar dynamical friction 
formula (Chandrasekhar 1942) 



dt 



X = 



-v M 

vm 
\/2<j v 



4tt In AG 2 (M + m)n m 



v 3 
"m 



erf(A) 



2X _x- 



(219) 



where a v is the velocity dispersion of the lighter particles. 
The Coulomb logarithm A arises from an integration over 
all impact parameters smaller than an upper limit Z max and 
is given by 



A ; 



J2, 



G(m + M) 



(220) 



Although encounters in the gravitational field of the sun de- 
viate from pure two-body scatterings, it is safe to neglect 
the presence of the sun if the encounter velocity is large 
compared to the Hill velocity 22 flRmn- Thus the classical 
dynamical friction formula is also applicable to planetesimal 
encounters in the high velocity regime, though a generalisa- 
tion to triaxial velocity distributions cti is necessary (Binney 
1977,see e.g.): 

dVM,i 

dt 



v M ,iV2iTG 2 \n(A) (M + m)n mBi (221) 



du 



(222) 



An additional complication is the choice of Z max (i. e. the 
choice of the Coulomb logarithm). There are several scale 
lengths, which could determine the largest impact parame- 
ter Z max : The scale height of the planetesimal disc, the radial 
excursion due to the excentric motion of the planetesimals 
and the Hill radius of the planetesimals. As it is not possible 
to derive a unique expression for Z max from first principles, a 
proper formula is often fitted to A?-body calculations (com- 
pare Eq. 235). The velocity dispersion of a planetesimal disc 
is triaxial with T^/T r = 1/4 and T z /T r « 1/4. We take 
these values and expand Eq. 221 for small velocities vm- 



dv 



M.r 



dt 

dvM,<t> 
dt 

dv M ,z 



dt 



-1.389 vm,. 



-3.306 »m. 
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V2^G 2 ln(A)(M + 


m)riom 



r.3/2 



(223) 



The derived expressions provide a compact tool to analyse 
dynamical friction in disc systems. However, the involved 
approximations are too severe compared to the needs of 
an accurate description. While these concise expressions are 
valuable for basic estimations, the following sections derive 
viscous stirring and dynamical friction formulae for a plan- 
etesimal system in a rigorous way. 



16.3 High Speed Encounters 

We return to the Boltzmann equation as a starting point for 
the derivation of the scattering coefficients: 



df 
dt 



+ v- V/- V$- 



dl 
dv 



(224) 



Eq. 10-12 provide the coordinate transformation. 



In virtue of the ansatz for the distribution function (see 
Eq. 210), it is sufficient to derive the time derivative of the 
second order velocity moments T r and T z . Since the distri- 
bution function is time independent in the absence of en- 
counters, only the collisional term contributes to the time 
derivative of the velocity dispersions Tk (k G (r, z, (j>) in the 



22 Whenever relative velocities are classified as "high" or "low" 
in the following sections, a comparison with the Hill velocity is 
implied. 
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following) : 



dpT k 
dt 



J d A vmv\ 



(225) 



coll 



The collisional term invokes the averaging over many differ- 
ent scattering trajectories and is, given that the underlying 
encounter model is analytically solvable, still too complex 
to derive an exact expression. If most of the encounters are 
weak - a realistic assumption in a planetesimal disc - it is 
possible to expand the collisional contribution in terms of 
the velocity change Av,. This is the Fokker-Planck approx- 
imation (see e.g. Binney 1994)) 



coll 



£_ [/D( a„,,] + 



It. 



2 ^ dvidv-j 



[D{A Vl ,A Vj )] 



(226) 



where the diffusion coefficients D contain all information on 
the underlying scattering process. Next we consider two in- 
teracting planetesimal populations m, m* with distribution 
functions 



2n 2 T r T z m 



exp 



v'i + nV 



r = 



27r 2 T*r z *', 



2T r 



2T Z 



2T* 



2T-* 



(227) 



to evaluate the terms in equation 226. We follow Stewart 
(2000) except some minor changes in the notation. The col- 
lisional term requires an averaging over the velocities of the 
two interacting planetesimals m and m* : 



d(pVk) 
dt 



= 2TTG 2 mm* 2 J d A v J d A v*ff*: 



2A(m + m*)u k v k Bu 2 + (2C - B)3u 2 k ) 



Wfe = Vk — v k 



A = ln(A 2 + 1) C= A2 - 1 B — A — C (228) 



■u-integration: 

d(pv 2 } 
dt 
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= 2nG 2 mm* 2 J d A w j d A uffx 
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All integrals are solvable and give the result 



G 2 p* 



d{Vr,4>) 

dt ~ 2 v / 2(T r + T r *) 3 / 2 

[B(r; + T r )m* JrAP) + 2A(T r *m* - mT r )H rA (P)} 

d(vl) = GV 
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[B(r: + T z )m*J z (P) + 2A{T* z m* - mT z )H z (/3)] 
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T r +T* 



(230) 



where six auxiliary functions are introduced to arrive at a 
more compact notation: 

a = y/l - 3s 2 6 = v/l - (1 - f3 2 )x 2 



f x 

H r := 8\/n / — dx 

Jo ab 

J a(/3a + b) 

H z := 80F / P ) ' dx 

J b(/3a + b) 



J<t> 
J z 



— —2H r + H^ + H z 

— H r — 2Ha, + H z 
= H r + H& — 2H Z 



(231) 



Since these are non-trivial functions, we apply a standard 
Chebyshev approximation for f3 € [0, 1]: 



1 

f(x) ^^c k T k {x) - -co 



(232) 



A coordinate transformation to the relative velocity u and 
the modified centre-of-mass velocity w 



w k = < 



V k + 



V k + 



(m*T* - mT r )u k 

(m + m*){T r +Tf) 
(m*T* - mT z )u k 



V = 



{m + m*)(T z +T;) 
mv + m*v* 
m + m* 



for k G {r, 0} 
for k — z 



(229) 



further simplifies the double integral. Thus the integration 
separates in a simple integral over to and a more demanding 



Table 9 summarises the Chebyshev coefficients. A final z- 
averaging yields the expressions: 



d{v 2 >(l> ) 



(233) 



dt 4^(T r + T r *)3/2(T Z + T;)V2 

[B(T; + T r )m*J rA {p) + 2A(T r *m* - mT r )H r>4> {!3)] 
« = G ' 2 ^* x f234) 

[B(T; + T z )m*J z (P) + 2A(T* z m* - mT z )H z {P)\ 

The determination of a proper Coulomb logarithm A leaves 
room for further optimisation. A careful comparison with 
A^-body models gives rise to the empirical choice (Ohtsuki 
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Jr 


J4, 


Jz 


H r 




H z 


co 


-10.34660733 


1.81674741 


8.52985992 


11.00434580 


6.94989422 


4.71219005 


Cl 


4.69990443 


2.95397208 


-7.65387651 


-2.64707927 


-2.06510182 


1.47084771 


C'2 


-1.25533220 


-1.18724874 


2.44258094 


0.60969641 


0.58700192 


-0.62294130 


C3 


0.30288875 


0.37775788 


-0.68064662 


-0.13815856 


-0.16311494 


0.18968657 


c 4 


-0.07040537 


-0.11070339 


0.18110876 


0.03112047 


0.04455314 


-0.05271757 


C5 


0.01540098 


0.02922947 


-0.04463045 


-0.00669979 


-0.01130929 


0.01331068 


A 


0.006 


0.015 


0.022 


0.0025 


0.0058 


0.0066 



Table 9. Chcbyshcv coefficients of the auxiliary functions and H^. 



et al. 2002): 



f2i?Hill 



(235) 
(236) 



Ohtsuki et al. 2002 also report a further improvement by 
setting B = A. 



16.4 Low Speed Encounters 

Encounters in the low velocity regime exhibit a wealth of dif- 
ferent orbits, as the solar gravity field perturbs the two-body 
scattering. Only a small subset of the trajectories represents 
simple, regular orbits like Tadpole or Horseshoe orbits 23 . 
Hence an examination of this velocity regime is done best 
with a numerical study of the parameter space by integrat- 
ing the equations of motions numerically (see Eq. 16). 

Ohtsuki et al. (2002) integrated a large set of planetesi- 
mal encounters and extracted fitting formula? that cover the 
low velocity regime. Their expressions for viscous stirring 
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The stirring rate of the radial velocity dispersion approaches 
a finite value for very low velocity dispersions, while the stir- 
ring rate for the vertical velocity dispersion drops to zero as 
the velocity dispersion decreases. This different behaviour of 
the two limits is due to the encounter geometry: If two plan- 
etesimals have zero inclination, they may still excite higher 
eccentricities during an encounter, but they remain confined 



23 The most famous example of such a regular orbit are the two 
saturnian moons Janus and Epimctheus which share nearly the 
same orbit. 



to the initial orbital plane preventing any excitation of in- 
clinations. 

The respective expressions for the dynamical friction 
rates are: 



dT r _ 

dt ~ 6(m + m*)(T r +T*) 
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dt ~ 6(m + m*)(T z + T*) 
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As the stirring rates are only valid in the low velocity regime, 
Ohtsuki et al. (2002) introduced special interpolation coef- 
ficients d. These coefficients tend to unity for very small 
velocity dispersions, and drop to zero in the high veloc- 
ity regime. Thus the interpolation formulae are properly 
"switched off" in the high velocity regime, so they do not 
interfere with the known high velocity stirring rates. 



16.5 Distant Encounters 

All formulas include only the stirring rates due to close en- 
counters, but non-crossing orbits also contribute to the over- 
all change of the velocity distribution. As these distant en- 
counters lead to small changes of the orbital elements, the 
problem is accessible to perturbation theory; see Hasegawa 
(1990) for a detailed treatment. Stewart (2000) integrated 
the perturbation solution over all impact parameters to de- 
rive the collective effect of all distant encounters: 
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(239) 



a accounts for the uncertainty in the smallest impact param- 
eter that is regarded as a distant encounter. While distant 
encounters are already included in the interpolation formula 
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of the low-velocity regime, we use the modified expression: 

(1 - Ci) 



dT r \ l, n ,2 (die 1 ) 



dt 



dist 



2(m + rn*) 2 



(JVs,diBt)(l - Ci) (240) 



Stewart (2000) omitted the change in the inclination, as it 
is small due to the encounter geometry. Nevertheless we de- 
rived the integrated stirring rate for completeness, which we 
give in 241 A close inspection of the integrated perturbation 
shows that the above formula is roughly a factor (i 2 ) + (i* 2 ) 
smaller than the corresponding changes in the eccentricity. 



16.6 Gas Damping 

The presence of a gaseous disc damps the velocity dispersion 
of the planetesimals and introduces a slow inward migration. 
Adachi et al. (1976) used the drag law Eq. 75 to approxi- 
mate 24 the average change of the orbital elements: 
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rjg is the dimensionless velocity lag of the sub-Keplerian 
rotating gaseous disc. 



16.7 Unified Expressions 

All expressions for the different velocity regimes are con- 
structed such that a smooth transition between the different 
regimes is assured. Thus, a simple addition of all contribu- 
tions yields already the unified expressions 
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which cover the full range of relative velocities. Although 
only two populations m and m* were assumed, Eq. 244 and 
Eq. 245 are readily generalised to a multi-mass system by 
adding a summation over all masses. 



16.8 Inhomogeneous Disc 

The preceding derivations assumed a homogeneous disc, 
which simplified the calculation, since the integration over 
all impact parameters needed no special precaution. A more 



sophisticated consequence is that the spatial density and the 
density in semimajor axis space are equal: 



E(r) = E(o) = So 



(246) 



Density inhomogeneities break this simple relation, as par- 
ticles at the same radial distance could have different semi- 
major axes, and particles with the same semimajor axis are 
located at different positions at a given time. While both 
representations are equivalent (i. e. describe the same sys- 
tem in different ways), we chose the density in semimajor 



axis space as the primary density 
derived as: 



The spatial density is 
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Likewise, T r and T z are also functions of the semimajor axis. 

Furthermore, an inhomogeneous surface density invali- 
dates the averaging over all impact parameters. Planetesimal 
encounter are most efficient for impact parameters smaller 
than a few Hill radii, so the derivation is still valid if the sur- 
face density is roughly constant on that length scale. How- 
ever, a planetesimal that is large enough will "feel" the spa- 
tial inhomogeneities or even generates density fluctuations. 
Hence it is essential to extend the validity of the averaged 
expressions to inhomogeneous systems. We use the averaged 
expressions 
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as a starting point (df r , z /dt excludes the surface density, as 
opposed to the averaged expressions). The (yet unknown) 
scattering contribution dT r , z /dt as a function of the impact 
parameter b is our starting point for a general expression for 
a varying surface density: 
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We restate Eq. 249 in terms of a weight function w(b) 
dT(a Q ) r ,z _ /dTr, z \ 1 
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L 

(251) 



The numerical solution of the Hill problem (see Eq. 16) 
gives some insight into how the weight function w(b) changes 
with the impact parameter. Fig. 12 illustrates the change in 
e 2 of the relative motion during an encounter of two plan- 
etesimals that were initially on circular orbits. While the 
details depend on the initial inclination and eccentricity as 
well as on the selected orbital element, all result share some 
basic features. Small (compared to the Hill radius) impact 
parameters allow for a horseshoe orbit and the change in 
the orbital elements is small except a change in the semima- 
jor axis. Intermediate impact parameters which lead to close 
encounters provide the strongest perturbation, but they are 
also more susceptible to complicated dynamics (compare the 



24 A formal expansion at e = 0, i = 0, r) g = is not possible, 
since the drag law involves the modulus of the relative velocity. 
Kary et al. (1993) corrected a missing factor 3/2 in Eq. 242. 



25 We denote S(a) also as "surface density" and refer to a as a 
radial coordinate. However, all formulae are precise in discrimi- 
nating both representations in r and a. 
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Figure 12. Change of the relative eccentricity e 2 due to an en- 
counter of two bodies initially on circular orbits. b/H is the impact 
parameter in units of the Hill radius. The plot was obtained by 
integrating Eq. 16. 



resonant structures in Fig. 12). As the gravitational attrac- 
tion drops with increasing distance, non-crossing orbits yield 
ever smaller perturbations with increasing impact parame- 
ter. Aside from this qualitative behaviour, it is very difficult 
to derive precise expressions. While the limit of high ve- 
locities is accessible through the two-body approximation, 
any general formula involves some empiric interpolation to 
cover the full parameter space (Rafikov 2003b, a,see the ap- 
proximations of). Therefore we decided to approximate the 
weight function such that the main features of the true 
weight function w(b) are reproduced. While this approach 
is less accurate, it provides better insight into the involved 
approximations. We expand the surface density under the 
integral in Eq. 250 and compare the expansion coefficients 
for w(b) and the approximation w(b) to derive constraints 
on the choice of w(b). The lowest non- vanishing order is: 
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/ can be interpreted as the width of the heating zone. Con- 
dition 252 inspires our choice of the weight function w(b) 
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where T^ 1 ' and are the radial velocity dispersions of 
the interacting radial bins and I is adjusted to the findings 
of Ida (1993). The advantage of the bell curve is that it has 
a discrete counterpart 
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which makes the weight function readily applicable to the 
summation on an equidistant radial grid with spacing Aa. 



16.9 Diffusion Coefficient 

We concentrated on the evolution of the velocity dispersion 
so far, but scatterings among planetesimals also change the 
semimajor axis of the disc particles, inducing a diffusive evo- 
lution of the surface density: 
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The diffusion coefficient D is related to the typical change 
in semimajor axis Aa and the timescale T^Body on which 
planetesimal encounters operate: 
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If we neglect the radial displacement during an encounter, 
the change in semimajor axis is solely due to the change of 
the velocity: 
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An average over all orientations of the velocity v and the 
velocity change Av yields the mean square change in semi- 
major axis: 
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This yields the mean diffusion coefficient 
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where the time derivatives of the velocity dispersions T r and 
T z are taken with respect to encounters. 



16.10 Coagulation Equation 

We already stated the coagulation equation for a multi-mass 
system: 



d d 
= ^mn(t,m) + —F m (t,rn) 



(260) 



F m = - J J ni(t,m 1 )n2(t,m 2 )a(m 1 ,m2)v lcl 

M Te d(m,mi,m2)dmidm2 (261) 

Since the vertical density profile of a planetesimal disc is 
specified by the known distribution function, we insert the 



32 P. Glaschke, P. Amaro-Seoane & R. Spurzem 



isothermal density profile (see Eq. 216) in the coagulation 
equation 260 and integrate over z: 



1 



(262) 
£i(f,mi) E 2 (i,m 2 ) 



^/2n(h(m 1 ) 2 + h(m 2 ) 2 ) mi m 2 
a(mi,m 2 )v le iM le d(m, mi,m 2 )dmidm 2 (263) 

E(m) is a short-hand notation for the differential surface 
density j^. Further integration over all masses gives the 
total mass balance: 
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The calculation of collisional cross sections is closely related 
to the underlying encounter dynamics. A homogenous sys- 
tem introduces no systematic perturbation, hence an en- 
counter is a pure two-body problem which is analytically 
solvable. Thus it is possible to derive the cross section with- 
out any approximation. Since encounters in the field of a 
central star deviate noticeably from the pure Kepler solu- 
tion, the cross sections are also modified. While the cross 
section in the high velocity regime reduces to the two-body 
formula (except minor corrections), the low velocity regime 
is explored best by numerical calculations. It is not appropri- 
ate to disentangle the different contributions in Eq. 263, but 
to combine the various terms to the collisional probability 
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(265) 



which can be easily deduced from the fraction of colliding 
orbits in Monte-Carlo simulations. An accurate expression 
for the collisional probability should include the two-body 
cross section in the limit of high velocities and the numerical 
data for the low velocity regime as well. We use numerical 
calculations from Greenberg (1991); Greenzweig (1992) 26 
as a basis for a unified fitting formula 
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which gives an effective cross section a for planetesimal- 
planetesimal encounters. Eq. 266 reduces to the well-known 
gravitational focusing formula in the limit of high velocities: 



(269) 
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(270) 



Ri + R2 

If the vertical velocity dispersion is small, the disc becomes 
two-dimensional and the cross section is proportional to R. 



26 Their work includes an averaging over the Raylcigh distributed 
inclinations and eccentricities of the colliding planetcsimals. 



The main differences to the two-body cross section 269 is a 
finite gravitational focusing factor, since the Keplerian shear 
inhibits a zero relative velocity, and a finite collisional prob- 
ability for very small velocities, again due to the shear which 
provides a finite influx of particles. 

The precise calculation of the coagulation kernel should 
include an integration over all semimajor axes with a proper 
weighting kernel. As collisions among particles in the statis- 
tical model play only a major role when the system is still 
homogenous, we omitted this contribution. In addition, this 
helps saving computational time, since the solution of the 
coagulation equation is very costly. However, interactions 
between JV-body particles and the statistical model include 
spatial inhomogeneities properly (see section 17). 

16.11 Collisional Damping 

Collisions are a dissipative process that removes kinetic en- 
ergy from the planetesimal system and damps the eccentric- 
ity. Low speed encounters leave the colliding bodies intact 
and damp the relative velocities through inelastic collisions. 
In contrast, high velocity encounters disrupt the colliding 
bodies and turn them into an expanding cloud of fragments. 
As a major part of the initial kinetic energy is converted 
into heat, the fragments disperse with rather low velocities 
thus reducing the overall velocity dispersion. We formulate 
the dissipation due to collisions analogue to Eq. 262: 
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Qred.fc is the kinetic energy redistribution function and Fq^ 
is the associated flux across the mass distribution, k indi- 
cates the two velocity dispersions. Q YC d,k is a complex func- 
tion, since the disruption of a planetesimal produces frag- 
ments with a large scatter in velocities and a complicated 
velocity field. The velocity of a fragment consists of two con- 
tributions: The ejection velocity relative to the target and 
the velocity v of the target within the corotating coordinate 
system. Owing to the strong dissipation, fragment velocities 
are dominated by the motion of the centre-of-mass of the 
two colliding bodies. Thus we neglect the ejection velocities 
and estimate the centre-of-mass motion. The initial kinetic 
energy of two colliding bodies is: 
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We separate E^i into centre-of-mass motion and relative 
motion and average over vi , v 2 : 
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Most of the relative kinetic energy is dissipated during the 
collision, so we neglect the relative motion after the incident. 
The final energy is 



(-Eflnai) _ m\T\ + mjTb 
(E ini ) ~ M(miTi + m 2 T 2 ) 



(274) 
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which gives the drift motion v of the expanding fragment 
cloud: 



V 2 = ^(Sflnal) 



(275) 



Qred is therefore coupled to the fragment redistribution func- 
tion M re< j 



Qrcd = V Af rc d + Qdis 



(276) 



where the additional function Qdiss removes the dissipated 
energy. 

16.12 Correlation 

The statistical model of a planetesimal disc does not only re- 
quire a large particle number to assure a proper description 
of the system by a distribution function, but also the un- 
corrected motion of the planetesimals. Each of the formulas 
derived before involves the averaging over different impact 
parameters to some extend, in combination with the vital 
assumption that all distances are equally probable. As long 
as all particles are subjected to perturbations by surround- 
ing bodies, strong correlations are suppressed. This applies 
to the early stages, but the formation of protoplanets intro- 
duces a few dominant bodies that are not susceptible to the 
perturbation of the field planetesimals. Orbit repulsion gives 
rise to a regular spacing of the protoplanets, which prevents 
mutual collisions. Therefore not all impact parameters are 
equally probable due to this strong correlation. Hence a sta- 
tistical model is inherently not applicable to the late stages 
of protoplanetary growth. 

Since statistical models are superior to TV-body calcula- 
tions with respect to speed and (effective) particle number, 
modifications have been proposed to remedy this problem. 

The statistical model by Wetherill (1993) uses the fol- 
lowing solution: A gravitational range Aa (or buffer zone) 
is attached to each planetesimal, which represents the min- 
imal spacing that allows for stable orbits. They propose the 
expression 



Aa = /a-Rhui + \/2T r /fP 



(277) 



where /a is the minimal spacing in terms of the Hill radius. 
The value of /a is adopted from Birn (1973), who derived 
the minimal spacing that allows for stable orbits: 



/a = 2%/3 



(278) 



Thus it is possible to define a minimum mass m sep by the as- 
sumption that all bodies larger than this critical mass main- 
tain a clear buffer zone: 
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/ is the area covered by the buffer zones (overlapping is not 
taken into account, therefore / > 1 is possible), normalised 
to the ring area. Planetesimals smaller than m sop can not en- 
force a minimum distance to their neighbours, as the whole 
disc surface is already covered by the buffer zones of the 
largest bodies. Owing to the regular spacing introduced by 
the buffer zones, planetesimals larger than the critical mass 
are not allowed to collide with each other. This approach 
has also been employed by Inaba et al. (2001), who adopted 
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Figure 13. Covered fraction fc as a function of / for simulation 
S1FB at T = 10 5 years. The protoplanets are already formed and 
grow oligarchically. 



/a = 10, which is the mean distance of protoplanets accord- 
ing to the orbit repulsion mechanism. 

To shed light on the proper gravitational range and the 
validity of this approach, we defined an additional quantity 
fc which is the true area (i. e. overlapping is handled prop- 
erly) covered by the buffer zones in terms of the total area: 

f°° dE 2na(Aa) c , 

fc = / —dm 

I „ am m 

" fJ1 sep 



fc < 1 



(280) 



Fig. 13 shows the covered fraction fc(m) as a function 
of the integrated buffer zones f(m) for one of our hybrid cal- 
culations. Both values for /a are included as well as the two 
limiting cases random placing and perfect ordering. Though 
we tested also other values of /a, a spacing of ten Hill radii 
proved to be the best choice. 

Our own experience with this method indicates that it 
works reasonably well and agrees with Inaba et al. (2001) 
who used the same technique. However, this modification in- 
cludes the regular spacing of the protoplanets in a prescribed 
way, so any exploration of later stages, like the initiation 
of orbit crossing, is not accessible through this approach. 
Therefore we use it for comparison purposes only, since the 
hybrid code (see section 17) provides a much more general 
framework. 



16.13 Discretisation 

All involved quantities are only functions of a and m. There- 
fore we introduce a two dimensional grid, where E, T r and 
T z are cell centred quantities. Fig. 14 summarises the defini- 
tion of the two-dimensional grid. Since the full planetesimal 
size range covers several orders of magnitude in mass, we 
chose a logarithmically equidistant discretisation in mass to 
cover the necessary mass range in a reliable way. The ra- 
dial spacing of the grid cells is equidistant. Thus the grid 
setup for the mass discretisation reads (TV grid cells from 
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Figure 14. Numerical grid. The arrows indicate transport of ki- 
netic energy (red), spatial transport of mass(green) and accretion 
(black). Non-neighbouring cells are coupled by the coagulation 
kernel and the radial interpolation kernel. 
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The grid spacing 5 controls the number of cells which are 
necessary to cover a specified mass range. As the evaluation 
of the coagulation equation scales with the third power of 
the number of grid cells, 8 should be as large as possible. If 
the flux integral (see Eq. 260) is approximated in a standard 
way 
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a spacing 5 much smaller than 2 is required to guarantee a 
sufficient accuracy 27 . However, it is possible to use a spacing 
of 2 if special precaution is taken. Spaute et al. (1991) ap- 
proximated the surface density with a power law, thus tak- 
ing the gradient with respect to mass into account. While 
they reached only a sufficient accuracy with further special 
adaptations, we use a more rigorous approach. A large spac- 
ing 5 reduces the accuracy, since the partial flux (Eq. 282) is 
strongly varying even inside one grid cell. Thus we rearrange 



this expression to identify the most important terms, as we 
can see in Eq. 283. 

The strongest varying contribution Fv is now approxi- 
mated by a power law with respect to mj\ 



Fv(m) « Fv{mj) 



(284) 
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Thus Eq. 284 is used to provide improved partial fluxes F>' 
and thus we come to Eq. 285. 

Since the fragment redistribution function is a piecewise 
power law, an analytical solution of the integral is possible. 
Eq. 285 gives reliable results even with a spacing 8 = 2. The 
time derivative of the surface density reads 



+ Fi 



(286) 



which assures the conservation of mass within numerical ac- 
curacy. 

16.14 Integrator 

All contributions to the evolution of the surface density E 
and the velocity dispersions T r and T z are summarised by 
the following set of differential equations: 
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The Laplace operator is approximated in accordance with 
the equidistant radial grid: 
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We chose the Heun method 28 as the basic integrator for the 
statistical model. It is a second order accurate predictor- 
corrector scheme (X is a vector containing all the above 
quantities): 

X p = X n + Atf(X n ) 

X n+1 =X P + ^At(f(X p ) - f(X n )) + 0(At s ) (289) 

The Heun method is readily extended to an iterate scheme, 
which is equivalent to the implicit expression: 

X n+1 =X n + ~At{f(X n+1 ) + f(X n )) (290) 

This adds stability to the method and allows the secure in- 
tegration of stiff configurations that may appear during the 



° The name of this method is not unique. Some texts denote it 
27 Ohtsuki et al. (1990) give a thorough analysis of the impor- as the modified Euler method. The Heun method is a special case 

tance of the resolution. of the Runge-Kutta methods. 
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Figure 15. Interplay between the ./V-body component and the 
statistical component of the hybrid code. Black arrows indicate 
mass transfer, red arrows exchange of kinetic energy and green 
arrows indicate spatial structuring, respectively. 



runaway accretion phase. In practice, three iterations are 
sufficient to guarantee a stable integration. As the diffu- 
sive part is discretised with a first order accurate formula 
(see Eq. 288) , the whole iterated scheme is equivalent to the 
Crank-Nicolsen method. We choose a global time step for 
the statistical model according to the expression 

Vm S cj,X € {E,T r ,T z }^ (291) 

where the hybrid code (see next section) applies an addi- 
tional discretisation in powers of two to achieve a better 
synchronisation with the iV-body code component. 



At = min 



is powerful enough to provide a complete and accurate de- 
scription of the planetesimal problem, since each method is 
confined to a certain range of the particle number. However, 
these restrictions are complementary in the sense that each 
method covers a regime where the other method fails. This 
intriguing relation stimulated the construction of a hybrid 
code which combines the benefits of both methods. 

The basic idea is to introduce a transition mass mtrans, 
which separates the two mass regimes. Particles with a lower 
mass are treated by the statistical model, whereas larger par- 
ticles belong to the iV-body model. Though both parts are 
clearly divided in different mass ranges, they are connected 
by various interdependencies: 

(i) Direct collisions between particles lead to a mass ex- 
change. One process is the accretion of small particles by 
./V-body particles, but agglomeration within the statistical 
model can also produce particles larger than the transition 
mass. This requires the generation of new iV-body parti- 
cles. Energetic impacts may erode larger particles, so a cor- 
responding particle removal is also required for consistency. 

(ii) Mutual scatterings among iV-body particles and 
smaller planetesimals transfer kinetic energy. While energy 
equipartition leads to a systematic heating of the smaller 
field planetesimals, a consistent treatment has to include 
both transfer directions. 

(iii) Accretion and scattering by the A-body particles in- 
duce spatial inhomogeneities or even gaps in the planetesi- 
mal component, if the particles have grown massive enough. 
Likewise, the small particles could induce some structure in 
the distribution of the iV-body particles. Since the spatial 
structure is dominated by the stirring from few protoplanets, 
we neglect the latter process. 

Fig. 15 summarises this brief overview of the interactions 
between the two code components in a schematic diagram. 
The following sections explain the implementation of each 
interaction term in more detail. 



17 BRINGING THE TWO SCHEMES 
TOGETHER: THE HYBRID CODE 

We introduced two different methods to solve the plan- 
etesimal growth problem. On the one hand, we modified 
NbodyGH — h, which has been used so far mainly for the sim- 
ulation of stellar clusters, to adapt it to the special require- 
ments of a long-term integration of planetesimal orbits. On 
the other hand, we developed a new statistical code with a 
consistent evolution of the velocity dispersion, the capabil- 
ity to treat spatial inhomogeneities and a thoroughly con- 
structed collision treatment. Neither of the two approaches 



17.1 Mass Transfer 

An iV-body particle accretes smaller particles in its vicinity. 
We already derived expressions which describe agglomera- 
tion within the statistical model, so it is manifest to apply 
these formulae to derive the accretion rate of an iV-body 
particle. 

Most of the material is accreted within the cross- 
sectional area a (see Eq. 267), but the finite eccentricity 
of an orbit extends the accessible radial feeding zone. Thus 
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we assign the following surface density to each particle 



E(o) = 



M 



2ira\/2iTl 
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I 2 = (t/tt+ ^a 2 e 2 +T r /Q 



2l 2 
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by smearing it out over its feeding zone. T r is the radial 
velocity dispersion of particles in the statistical model with 
semimajor axis a. This density distribution is projected onto 
the radial grid to calculate the accretion rate. As the time 
step of the statistical model is much larger than the regu- 
lar step of an iV-body particle, the particle mass update is 
synchronised with the statistical integration. The projection 
technique allows the calculation of the accretion rates in a 
simple way, which gives the right size of the feeding zone 
and the proper total accretion rate. 

Particle generation is included in the following way: A 
"virtual" mass bin is introduced as the boundary between 
the statistical grid (denoted by the dashed area in Fig. 15) 
and the iV-body component. Its sole task is to store mass 
and kinetic energy that drives the statistical model towards 
higher masses. If the mass content exceeds one m trans , a new 
particle is created with inclination and eccentricity accord- 
ing to the stored velocity dispersions. 

The masses of the iV-body particles are regularly 
checked to detect any particle which dropped below the tran- 
sition mass. While this procedure would remove the particle 
and transfer the associated quantities back to the grid, we 
never observed such a particle erosion. 



17.2 Disc Excitation 

The projection of an iV-body particle onto the grid with 
the help of a proper weight function is also useful for the 
calculation of the disc excitation due to stirring by the larger 
particles. Since the Hill radius sets the proper length scale 
for planetesimal encounters, the weight function is modified 
to 
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where T r is the radial velocity dispersion of the heated plan- 
etesimal component. The velocity dispersion of the stirring 
iV-body particle is (in accordance with Eq. 11 and Eq. 12): 



f r,0 



^(na ) 2 e 2 T z , 



-(S2oo) I 
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We employ the orbital elements as mediators between the 
fast varying instantaneous position and velocity of a parti- 
cle and the slow evolution of the statistical model, which 
operates on a longer relaxation timescale. In virtue of the 
projection of the particle, we readily apply the standard in- 
teraction terms (see section 16) to evaluate the additional 
heating due to the presence of iV-body particles. 

17.3 Pseudo-Force 

While an iV-body particle is moving through the disc, it also 
interacts gravitationally with the particles in the statistical 
model. The collective effect of all these encounters leads to 



a change in the orbital elements of the JV-body particle. 
Again, we project the iV-body particle onto the grid and 
evaluate the stirring rates T r and T z , which correspond to a 
change in the orbital elements: 
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These time derivatives of eccentricity and inclination are 
translated to a pseudo-force that effects the desired change 
of the orbital elements. We chose the ansatz 



■ x.y — Cr(y x ,y (Vie)x,y) 



C z v z 
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where vk is the local Keplerian velocity. In addition, wc 
tried a simpler expression 

r ■ v 



Fx y — 2C r T x y ; 

F, = C z v z 



(298) 



without any significant differences in the accuracy or the 
simulation outcome. The proper friction coefficients are: 



Cr — 

c z = 



2T r 
2T Z 



(299) 



Since the relevant quantities are the time derivatives of the 
orbital elements, any other pseudo-force is also applicable. 
Though this approach yields the right mean change of the 
orbital elements, it lacks the statistical fluctuations from 
the particle disc. Hence the distribution of the orbital ele- 
ments of the iV-body particles is artificially narrowed, which 
is especially important when the iV-body particles and the 
statistical particles have a comparable mass. As the mass 
contrast between the two code parts is quite significant in 
planet formation simulations, it is safe to neglect the fluctu- 
ating part without major restrictions on the realism of the 
simulations. 

The friction coefficients d are kept constant between 
two integration steps of the statistical model. While a more 
frequent update of the coefficients would be easily possible, a 
regular update on the basis of the statistical time step is ac- 
curate enough. Moreover, each update poses a considerable 
computational effort (roughly equivalent to 1000 force eval- 
uations), so our approach also saves valuable computational 
time. 



17.4 Spatial Structure 

The first insight into planetesimal formation was obtained 
by the particle-in-a-box method, which invokes the under- 
lying assumption that the planetesimal disc stays homoge- 
neous throughout the protoplanet growth (Greenberg et al. 
1978, see e.g.). While few large bodies introduce some coarse- 
graininess of the surface density, all smaller bodies are as- 
sumed to be evenly spread in the disc. Research on the in- 
teraction of protoplanets showed that this is an oversimpli- 
fication, as bodies that are massive enough could open gaps 
in their vicinity (Lin & Papaloizou 1979; Rafikov 2001, see 
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Table 10. Parameters of the statistical and the TV-body gap simulation. The perturber is placed at the centre of the ring. 
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Figure 16. Gap opening in a planetesimal disc. The gap is fully 
developed after 2000 years. Table 10 summarises the initial con- 
ditions for the comparative runs. 



Figure 17. Mean square eccentricity and inclination of the 
smaller planetcsimals in terms of the reduced Hill radius H of 
the protoplanet according to simulation Gl (TV— body) and G2 
(Statistic). 



e.g.). Gap formation does not only change the overall surface 
density, but also controls the accretion onto the protoplanet 
through the amount of planetesimals in the feeding zone. If 
gap formation is too effective, the growth of the protoplanet 
may well stop before the isolation mass is reached. Hence 
any hybrid code should provide a framework that allows 
this mechanism to operate. A necessary condition is a radial 
density grid with a sufficient resolution to describe possibly 
emerging gaps. A too low resolution suppresses local pertur- 
bations from the protoplanets by a simple averaging, thus 
inhibiting the formation of any spatial inhomogeneities. A 
second requirement is that the interaction terms relating 
statistical model and TV-body model include the local inter- 
action between particles and the statistical component in a 
proper way. 

Our hybrid approach includes gap formation implicitly 
through the diffusive terms. A protoplanet heats only the 
planetesimals in its vicinity (defined by the heating kernel), 
thus also increasing locally the diffusion coefficient. Hence 
the surface density drops due to outward diffusion of the 
planetesimals, given that the protoplanet is massive enough. 
The minimum gap opening mass is set by the condition that 
the protoplanet controls the random velocities of the field 
planetesimals in its heating zone (see Eq. 30), which is equiv- 
alent to the independently derived gap formation criterion 
(compare Eq. 32). 

Although our algorithm invokes a simplified picture of 
the protoplanet-planetesimal interaction, it is surprisingly 
accurate with respect to the width of the forming gap and 
the opening criterion. Fig. 16 shows a simulation which ex- 



amines the accuracy of our approach. The overall perfor- 
mance of the statistical code is quite remarkable, except a 
significant overestimation of the surface density at the gap 
boundary compared to the TV-body model. This deviation 
is due to the improper treatment of strong planetesimal- 
protoplanet encounters, which exceed the diffusive approxi- 
mation. Moreover, the higher concentration of planetesimals 
near the gap boundary leads to an additional overestimation 
of the velocity dispersion of the smaller planetesimals in the 
statistical calculation (see Fig. 17). While the comparison 
with the TV-body calculation clearly indicates a necessary 
improvement of the treatment of spatial inhomogeneities, 
our approach catches the main features of gap formation. 



IT. 5 Transition Mass 

Since the inventory of the new hybrid code is now completed, 
we turn to the specification of the transition mass mtrans- 
The mass boundary between statistical and TV-body part 
has a major influence on the realism and the speed of the 
simulation. On the one hand, optimisation with respect to 
speed favours a large transition mass, whereas a reasonable 
resolution of the transition between the two components in- 
troduces some upper limit. 

Hence we identify first the set of large masses, which 
controls the velocity dispersion of the disc, since these ob- 
jects are also possible candidates for gap opening. The in- 
spection of all involved stirring terms gives approximately 
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the inequality: 
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While this is a necessary condition to select all potential ma- 
jor perturbers, criterion 300 does not imply that all particles 
in the selected mass range exert indeed a strong influence 
on the disc. The number of possible gaps - and therefore the 
number of perturbers associated with them - is ultimately 
limited by the available space. Thus we integrate the area 
of all potential gaps (width « /a^hhi) and normalise it to 
the total disc area: 



fc 



JA-, dm (301) 

dm m 



If the covered fraction fc is much larger than one, it is 
possible to increase the transition mass until the condition 



(302) 



is fulfilled. Of course, condition 300 and 302 defined only an 
upper limit of the transition mass, so the adaptation of a 
lower value is also possible. Though there are two reliable 
conditions at hand, the transition mass is still a function 
of time owing to the time evolution of the density S(m). 
Therefore we chose a priori a fiducial value of the transition 
mass, run the simulation and conduct an a posteriori check, 
whether the initial choice matches our requirements at any 
evolutionary stage of the disc. A reliable value for a solar 
system analogue at 1 AU is 



TOtrans ~ 3 X 10 "Mq 



(303) 



which restricts the number of iV-body particles to a 
tractable amount. Later stages would allow an even larger 
transition mass, but the current hybrid code does not in- 
clude any dynamical adjustment of the transition mass at 
runtime. 



17.6 Boundary Conditions 

Any numerical simulation is limited to a finite simulation 
volume and a finite time interval. Therefore it is mandatory 
to introduce proper boundary conditions which provide a 
reasonable closure of the simulation volume. 

While boundary conditions with respect to time are the 
familiar initial conditions, the choice of the spatial bound- 
ary conditions for the various involved quantities depends on 
the problem at hand and the type of the boundary. A sim- 
ulation boundary can be due to physical reasons (like walls 
of a concert hall, surface of a terrestrial planet) or simply 
due to a limitation in computational power that inhibits the 
complete numerical coverage of the problem. 

The current capability of the hybrid code sets limits on 
the radial range as well as on the covered mass range, which 
a simulation can handle in a reasonable time. Hence we have 
to introduce artificial boundaries in radius, and a lower limit 
for the mass grid. 

Any migration process couples the evolution of a lo- 
cal ring area in the planetesimal disc to the evolution of 
the whole disc. Inward (or outward) migrating material also 
transports information on the radial zone where the material 
originated from. As this information is not available within 



the frame of a local simulation, any choice of the boundary 
condition alters the evolution to some extend. 

However, we focus on a formation stage where migra- 
tion is not a dominant process, but provides only removal 
of the smaller collisional fragments. Thus we apply closed 
boundary conditions for the outer and inner radius of the 
ring area (i.e. all fluxes vanish at the boundary), and an 
open boundary for the lower end of the mass range. While 
these conditions exclude the study of migrational processes, 
we gain clearer insight into the protoplanet growth. 



18 DISCUSSION AND CONCLUSIONS 

The formation of planetary systems represents a challenge 
from a numerical standpoint. The dynamical problem spans 
over many orders of magnitudes in length and demands the 
combination of different techniques. We have presented a 
composite algorithm that brings together the advantages 
of direct-summation tools and statistics for the description 
of the planetesimal disc. Direct-summation A— body tech- 
niques have been around for some decades and have proven 
their accuracy in a very large number of studies of stellar 
clusters such as galactic nuclei and globular and open clus- 
ters. We deem it to be the numerical tool to integrate the 
motion of the bodies for the very precise integration of the 
orbits and treatment of close encounters. Typically, in a sim- 
ulation of a stellar system, the energy is conserved in each 



timestep by E/AE 



10" 



(where E is the total energy 



and AE the difference between the former and current to- 
tal energy for a specific time), so that even if we integrate 
for a long time the cluster, the accumulated energy error 
is negligible. Nevertheless, porting the numerical tool to the 
problem of planetary dynamics is not straighforward and re- 
quires important modifications and additions. In this work 
we present them in detail: the neighbour radius selection for 
the protoplanets, the Hermite iteration and we introduce 
for the very first time the new extended Hermite scheme, 
since the usual Hermite scheme is not sufficient to integrate 
planetesimal orbits accurately enough. Then we bring in new 
forces to the problem, namely the introduction of the central 
potential of the star, as well as the drag forces, which depend 
on the gas density and size of the planetesimals. Hence, the 
regularisation scheme, crucial to exactly integrate the close 
encounters, has to be accordingly modified. We then intro- 
duce the disc geometry and discuss the required changes to 
the neighbour scheme and prediction, as well as the commu- 
nication algorithm and block size distribution. 

For the statistical description of the planetesimal disc 
we employ a Fokker-Planck approach. We include dynamical 
friction, high- and low-speed encounters, the role of distant 
encounters as well as gas and collisional damping and then 
generalise the model to inhomogenous discs. We then de- 
scribe the combination of the two techniques to address the 
wole problem of planetesimal dynamics in a realistic way via 
a transition mass to integrate the evolution of the particles 
according to their masses. 

In particular, we introduce and describe the extended 
Hermite scheme, which reduces the the energy error by three 
orders of magnitude with the same number of force evalua- 
tions, compared to the standard version of Nbody6+ + . 

While the implementation and some code details are 
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Table 11. Parameters of all test simulations. The transition mass in T4b is mtrans = 3.1 X 10~ 10 Only simulations T3, T4a-T4c and 
T5 include collisions. All values are scaled to M c = G = ro = 1. 



newly introduced to the field of planet formation simula- 
tions, the first hybrid approach was developed in the early 
90's. Spaute et al. (1991) (further improved in Weiden- 
schilling et al. 1997) constructed a hybrid code with a statis- 
tical component to treat the smaller particles and a special 
treatment for the larger particles. A statistical model cov- 
ers the field planetesimals with the help of a distribution 
function (similar to Wetherill (1989)), whereas the larger 
particles are individually stored and characterised by mass, 
semimajor axis, eccentricity and inclination. While the inter- 
action between these single particles and the statistical com- 
ponent is expressed by standard viscous stirring and dynam- 
ical friction terms, perturbations among the single particles 
are equated in a different way. First, the probability of an 
encounter of two neighbouring particles is calculated. This 
probability is used in a second step to decide whether a (nu- 
merically integrated) two-body encounter of the neighbour- 
ing particles is carried out to derive the change in the orbital 
elements. Though these two well-defined code components 
justify to speak about a hybrid approach, the Monte-Carlo 
like integration of the largest particles is still closely related 
to a statistical treatment. 

A modified A^-body approach is used in the work of Lev- 
ison & Morbidelli (2007). Their method covers the largest 
particles by a direct A^-body code, which includes the 
smaller particles as "tracer" particles. The term "tracer" 
indicates that each particle represents a whole ensemble of 
planetesimals. In a similar line of approach and inspired by 
this idea, Levison et al. (2010) modified a symplectic algo- 
rithm, Symba, to study the formation of giant planet cores. 
However, they made some assumptions in order to calcu- 
late the gravitational interaction between the planetesimals. 
In particular, they ignored totally close encounters between 
planetesimals. 

Ormel & Spaans (2008) present in their work a scheme 
based on Monte Carlo techniques to cover the vast range of 



sizes. For this, they assign more resolution to those parti- 
cles that are more relevant to the interactions, typically the 
largest bodies. Smaller particles are grouped and treated 
collectively, which means that they all share the same mass 
and structural parameters. This classification is done in ac- 
cordance to the "zoom factor" , a free paramenter. Later, 
Ormel et al. (2010a) presented an detailed comparison of 
their Monte Carlo code with other techniques, in particular 
with pure direct-summation N— body results and other sta- 
tistical studies and found that system leaves the runaway at 
a larger radius, in particular at the outer disc. With their 
simulations, the authors propose a new criterion for the run- 
away growth-oligarchy transition: from several hundreds of 
km in the inner disk regions up to a thousand km for the 
outer disc (Ormel et al. 2010b). 

Bromley & Kenyon (2006) published a description of a 
hybrid method with a basic approach similar to our work. 
They employ two velocity dispersions and the surface den- 
sity of the planetesimals to describe the planetesimal system. 
The statistical component includes migration of the plan- 
etesimals and dust particles due to gas drag and Pointing- 
Robertson drag. In contrast to our approach, they did not 
include mass transport due to the diffusion of the planetesi- 
mals, which precludes the study of spatial structures induced 
by the protoplanets. One must note also that their method 
uses the standard discretisation of the collisional flux (see 
Eq. 282) and thus restrict the spacing factor to 8 < 1.25 
(Kenyon 1998). Bromley & Kenyon (2006) chose a set of 
test calculations which focused less on the technical aspects 
of their method, but on an overall comparison with a se- 
lected set of standard works on planet formation. Their test 
simulations are in good agreement with the references simu- 
lations, thus indicating a comparable quality of the method. 
Four years later, the authors presented an updated version 
of their code for planet formation. The new characteristics 
of the code included ID evolution of the viscous disc, gas ac- 
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cretion on to massive cores, as well as accretion of small par- 
ticles in planetary atmospheres (Bromley & Kenyon 2010). 

While a variety of hybrid approaches emerged over the 
past years, this technique is still far from a routinely appli- 
cation and is still challenged by many open issues. Hybrid 
codes bear the potential to address the dynamical evolution 
of a whole planetary system, the later stages of protoplanet 
formation initiate a strong interaction with the gaseous disc, 
which may require more diligence than the inclusion of a few 
additional interaction terms. However, the development is 
picking up speed, which places our work in a good position 
for further research. 



APPENDIX A: CENTRAL FORCE - 
DERIVATIVES 

Central force F per mass (i. e. acceleration) and its time 
derivatives are: 
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The F^ denote the central force and its time deriva- 
tives, whereas a and a refer to the total acceleration of the 
particle. The assumption that x, v, a and a are independent 
of each other allows the derivation of averaged expressions 
for particle-particle interactions: 
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We combine these expressions with Aarseth's time step for- 
mula to derive the regular time step as a function of the 
neighbour sphere radius R s : 
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r is the average particle distance and r c i oso is the impact 
parameter for a 90-degree deflection. 



APPENDIX B: SCALABLE COLLISIONS FLUX 

The mass flux according to the perturbation equation 161 
is: 

F p — — J j (n(m2)An(mi) +n(mi)An(mj)) 
a{m{)v{m2)m\f m {m\/m 1 e)dm 1 dm,2 

= F W + F (2) 

Firstly, we employ the substitution 
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to solve for the partial flux F^- 1 ' 
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The second contribution F^ 2 ' requires a slightly different 
transformation: 



-V(i+") 



mi = mx\e 

( mil \ T 



i ■ .. 

+ 2„ 



m 2 — mo 



\ m J 
Thus the partial flux F (2) is: 
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We change to a new set of logarithmic coordinates 
u = ln(m/mo) ui = ln(xi) § = 

l + a 

which transforms the total flux F m to a convolution integral: 
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F p = -nlmlaovo J [g(u + u 1 )Gi(u 1 )+ 

gipiu + m +si))G 2 (ui)]diii (B7) 

< B8 > 

p = 1 refers to the already derived solution for self-similar 
collisions. Hence we expand Eq. B7 at p — 1 and retain only 
the zeroth-order moment of the fragmentation kernel: 

F p = -n 2 mocr « [#(m)Gi,o + (fl(u)+ 

„(p-l)^)G 2 ,o] (B9) 
This expression is equivalent to 

F p = —nlmlaoVo(g(u)Gifi+ 

[g(u) + (p- l)(g(u) - <?(0))]G 2 , ) (BIO) 

where higher derivatives of g(u) are neglected. Hence we 
recover the same functional form of the perturbed mass flux 
F v as for self-similar collisions: 



Fp = —n m croVog(u) (Gi,o + pG 2 ,o) + const. 
ocS~ k ' 



(Bll) 
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APPENDIX C: COAGULATION EQUATION 

While the success of a general approximation of the coag- 
ulation equation depends heavily on the used coagulation 
kernel, we nevertheless provide a more general approach to 
embed section 10 in a broader context. The standard coag- 
ulation equation is: 

d d 
= —-mn(t,m) + ■ 7 —F m (t,rn) 
at dm 

F m = — J I n(t,m 1 )n(t,m 2 )cT{m 1 ,m2) 

u re iM red (m, mi,m2)dmidm,2 (CI) 

In virtue of our experience drawn from the perturbation 
expansion, we transform the coagulation equation to loga- 
rithmic coordinates 



u — ln(m) 



(C2) 



and employ the size distribution g(u) relative to the steady- 
state solution n eq (m): 

= ^g(u,t)n cq (u)e 2u + —F u (t,m) 

F u = — J j g{t,ui)g{t,U2)K(u,ui,U2)duidu2 (C3) 

K(u, U\,U2) is the properly transformed new coagulation 
kernel. g(u) is expanded under the integral to arrive at a 
moment expansion of the flux F u : 



F u = -K 00 {u)g(uf - {K 10 (u) + K 01 (u))g{u)^ + . 



Ki. 



J J K(u,m, 



U2)U\u\du\du2 



(C4) 



Retaining only the leading order terms, we recover an ap- 
proximate coagulation equation which is similar to the in- 
viscid Burgers' Equation 29 : 

= jLg(u, t)n oq (u)e 2 " - |- (K 00 (u)g(u) 2 ) (C5) 
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